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SUMMARY 

. . 

|An  analytical  model  is  presented  for  predicting  the  transient  one- 
dimensional  thermal  performance  pf  a  charring- ablator  heat-protection  system 
when  exposed  to  a  hyperthermal  environment.  The  heat-protection  system  is 
considered  to  consist  of  an  ablation  material  and  backup  structure.  The 
ablating  material  is  further  considered  to  consist  of  three  distinct  regions 
or  zones:  char,  reacting,  and  virgin  material. 

A  FORTRAN  IV  digital  computer  program  (STAB  II)  utilizing  an  implicit 
finite  difference  formulation  has  been  written  for  the  IBM  "JO^b/kO  computer 
system.  The  program  considers  one  ablating  material  and  a  maximum  of  12  back¬ 
up  materials  with  conduction  or  radiation  and/or  convection  allowed  between 
materials.  Thermal  properties  of  all  materials  are  temperature  dependent, 
with  the,  properties  of  the  charring  material  also  being  state  dependent. 

The  governing  differential  equations  and  their  implicit  finite  difference 
formulation  are  presented.  The  program  input  and  output  are  described  in  detail. 
The  FORTRAN  program  statements  and  nomenclature  are  presented.  Also,  the 
theoretical  and  experimental  results  are  compared^ 

INTRODUCTION 


^The  analysis  and  design  of  thermal  protection  systems  for  entry  into  an 
atmospheric  environment  have  resulted  in  a  voluminous  amount  of  literature  on 
the  general  s^ject^of  ablation^A^CSee  refs.  1  and  2  for  a  survey  of  information 
on  ablation.)  ’’The  [ablation  materials  may  generally  be  classified  into  three 
categories:  subliming,  melting  and  vaporizing,  and  charring.  The  charring 
ablator  normally  provides  the  most  efficient  thermal-protection  shield  for  the 
major  portion  of  a  manned  entry  vehicle.  This  report  describes  a  method  for 
predicting  the  thermal  response  of  a  typical  charring-ablation  material.  The 
response  of  a  charring  material  to  a  hyperthermal  environment  is  extremely 
complex,  and  the  mathematical  model  presented  to  analyze  the  transient  behavior 
of  the  material  contains  simplifying  assumptions  and  approximations  necessary 
to  afford  even  a  numerical  solution. 

^ /24'J — 


The  equations  derived  in  this  analysis  have  been  programmed  in  FORTRAN 
IV  for  an  IBM  'JOSh/hO  computer  system.  The  numerical  formulation  of  this 
digital  program,  designated  STAB  II,  is  such  that  an  implicit  solution  is 
obtained.  The  thermal  response  of  a  typical  charring  material  as  predicted  by 
STAB  II  is  compared  with  arc  tunnel  results. 

A  sample  problem  is  presented  in  appendix  A.  Program  usage  instructions, 
including  definitions  of  the  input  terminology,  are  presented  in  appendix  B. 
Appendixes  C  and  D  are  the  program  FORTRAN  IV  statements  and  definitions  of 
the  program  terminology.  A  general  flow  chart  of  the  program  is  presented  in 
appendix  E. 
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SYMBOLS 


collision  frequency 
specific  heat 

activation  energy 
exterior  view  factor 

view  factor— emissivity  product  to  cabin  environment 

heat  of  virgin  material  degradation 

total  enthalpy 

wall  enthalpy 

enthalpy  of  air  at  300°  K 

film  coefficient  between  backup  materials 

film  coefficient  between  last  backup  material  and  cabin 
environment 

thermal  conductivity 

mass  loss  rate  of  char  material 

gas  ablation  rate 

number  of  nodes  in  ablation  material 
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A9 
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order  of  reaction 

net  heat  rate  into  front  surface 

hot  wall  convective  heat  flux  with  blowing 

heat  flux  due  to  combustion 

cold  wall  convective  heat  flux  without  blowing 

radiation  heat  flux 

universal  gas  constant 
surface  recession  depth 
surface  recession  rate 

temperature  of  node  at  beginning  of  time  step 

cabin  environment  temperature 

temperature  of  node  at  end  of  time  step 

radiation  heat  sink  temperature 

thickness  of  ablation  material 

distance  from  surface  to  any  point 

heat  of  combustion  per  unit  weight  of  char 

thickness  of  a  node 

time  step  (0'  -  9) 

emissivity  of  material 

transpiration  cooling  efficiency 

initial  time 

final  time 
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transform  for  the  ablation  material 


p  density 

a  Stefan- Boltzmann  constant 

i|r  blocking  effectiveness  function 


Subscripts: 

c  charred  state 

i  node  number 

j  material  number 

v  virgin  state 


PROGRAM  DESCRIPTION 


The  following  general  requirements  were  established  before  writing  a 
digital  computer  program  to  analyze  a  charring  ablation  system: 

(1)  Stability  of  the  equations  for  all  applications. 

(2)  Machine  running  time  short  enough  to  make  use  of  the  program 
economically  feasible  (a  minimum  of  turn  around  time  per  problem). 

(3)  A  minimum  of  input  per  problem. 

(4)  A  wide  variety  of  boundary  conditions  for  application  to  both 
trajectory  data  and  ground  or  flight  test  data  analysis. 

STAB  II  has  been  formulated  in  FORTRAN  IV  to  analyze  the  transient 
thermal  performance  of  a  charring  ablator  heat  protection  system.  The  program 
considers  one  ablating  material  and  up  to  12  different  backup  materials  with 
or  without  air  gaps.  Pure  conduction  or  radiation  and/or  convection  between 
backup  materials  is  allowed.  The  ablation  material  may  be  divided  into  a 
maximum  of  50  nodes,  and  each  backup  material  may  be  subdivided  into  a  maximum 
of  10  nodes.  The  thermal  properties  of  the  materials  are  in  tabular  form  and 
are  temperature  dependent.  The  ablation  material  is  also  dependent  upon  its 
state,  that  is,  fully  charred,  partially  charred,  et  cetera. 

The  following  surface  boundary  condition  options  are  provided: 

(l)  Cold-wall  convective  and  radiative  heat  flux  tables  as  a  function  of 
time.  These  components  are  specified  separately,  since  mass  transfer  at  the 
surface  blocks  part  of  the  convective  heating  but,  in  general,  has  no  effect  on 
the  radiant  heating. 


(2)  Surface  temperature  as  a  function  of  time. 

(3)  Surface  recession  as  a .function  of  temperature  or  time.  Surface 
recession  as  a  function  of  temperature  and  pressure  is  also  available. 

Heat  loss  to  the  interior  environment  for  the  last  node  of  the  backup 
structure  can  be  specified  by  two  methods: 

(1)  Conduction  into  the  node  and  radiation  and/or  convection  loss  to  the 
interior  environment. 

(2)  Conduction  into  the  node  and  adiabatic  wall. 

The  STAB  II  numerical  formulation  of  the  equations  describing  the  response 
of  the  heat  shield  is  such  that  an  implicit  solution  has  been  obtained.  It  is 
well  known  that  numerical  solutions  of  partial  differential  equations  are 
subject  to  several  different  types  of  errors.  The  first  of  these  is  the  trun-» 
cation  error,  due  to  the  use  of  a  finite  subdivision.  This  error  maybe  reduced 
by  simply  choosing  a  smaller  subdivision,  AX.  The  exact  values  are  approached 
more  and  more  closely  as  AX  decreases.  The  second  kind  of  error  is  the  nu¬ 
merical,  or  roundoff  error.  The  way  in  which  this  numerical  error  grows  or 
decays  with  time  determines  the  stability  of  the  difference  equations. 

To  illustrate  the  differences  in  the  explicit  and  implicit  equation  form, 
consider  a  nonablating  homogeneous  solid.  The  one- dimensional  Fourier  conduc¬ 
tion  equation,  neglecting  any  heat  generation  terms,  is 


The  finite  difference  form  of  equation  (l)  written  in  the  conventional  forward 
time  step  or  explicit  form  for  the  i  node  is 


(Ti-i  -  Ti)  (^rfi±i) 


AX 

2k. 


T. 

x_ 

AX 

2k. 


AX 


2k 


1+1 


(2) 


where  the  prime  superscript  denotes  values  at  the  end  of  the  time  step 

A6  =  9'  -  0 

For  explicit  conduction  solutions,  the  following  stability  criterion  has 
been  established: 

PCP  Mf  >  2 

k  A9  ~ 
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which  places  an  upper  limit  on  the  time  step  A9  for  a  fixed  truncation  error. 
This  criterion  can  require  a  prohibitive  amount  of  machine  time. 

Liehmann  (ref.  3)  advocated  a  solution  of  the  equation  which  does  not 
require  this  stability  criterion.  The  finite  difference  equations  are  written 
in  a  backward  time  step  form  which  affords  an  implicit  solution. 

The  implicit  (backward  time  step)  difference  form  of  equation  (l)  for  the 
i^  node  is: 


Equation  (3)  uses  the  temperature  differences  at  the  end  of  the  finite  time 
interval  instead  of  the  beginning,  as  in  the  explicit  method  of  equation  (2). 
The  only  known  temperature  in  equation  (3)  is  T .,  but  there  are  corresponding 

equations  for  each  point  in  the  system,  and  all  are  solved  simultaneously  to 
yield  the  temperature  at  each  node. 

Collecting  all  unknown  temperatures  on  the  left  side  of  the  equation  and 
the  known  temperature  on  the  right  side,  equation  (3)  becomes 


Equation  (4)  is  of  the  form 

AT!  n  +  BT|  +  CT|  ,  =  D  (5) 

l-l  l  l+l 

STAB  II  generates  such  an  equation  for  each  node  in  the  system. 

Since  radiation  is  an  important  mode  of  heat  transfer  in  charring  ablative 
systems,  a  problem  is  encountered  in  any  equation  containing  a  radiation  term. 
The  radiation  heat  flux,  written  in  a  backward  difference  form  is: 

■  F6°  -  T-1‘)  (O 
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This  term  cannot  he  used  in  an  implicit  solution  since  the  unknown  temperature 
T^  is  to  he  the  4th  power.  The  4th  power  unknown  can  he  eliminated  hy  the 

following  linearizations: 


where 


AT  =  T'  -  T 
i  i 


let 


and  rewrite  equation  (7)  as 

(Ti)4  =(Ti)4(l  +  z)4 

If  Z  has  an  absolute  value  near  zero,  the  following  is  true 

(l  +  Z)^  S  1  +  4z 


Now  substituting  equation  (9)  into  equation  (8) 

(Ti)4  =(Ti)V  +  te)  + 


~  4T  ^T '  -  5T  ^ 
i  i  J  i 


(7) 


(8) 

(9) 


(10) 


Equation  (10)  is  a  linearized  approximation  of  equation  (7)  in  which  the  unknown 

temperature  is  only  to  the  first  power.  The  assumption  in  equation  (10)  is  that 

AT/T,  has  an  absolute  value  near  zero.  Figure  1  is  a  plot  of  the  error  obtained 

1  4  * 

when  (l  +  4z)  is  substituted  for  (l  +  Z)  .  For  most  ablation  problems  in  which 

the  surface  temperature  is  high  and  the  radiation  losses  are  significant,  the 

value  of  AT/T^  can  easily  be  controlled  to  values  of  less  than  ±0.1. 


Therefore,  equation  (6)  can  now  be  written 


Srad 


=  F6ff|4T  -TP’ 


(11) 
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Using  the  linearized  approximation  for  the  radiation  terms,  the  resulting  system 
of  implicit  difference  equations  constitute  a  tridiagonal  matrix  of  the  follow¬ 
ing  form: 


B,  Tn 

+  C,  T„ 

1  1 

1  2 

A2T1 

+  B2T2 

+  C2T3 

A,T_ 

+  B,T, 

3  2 

3  3 

D, 


= 


=  D, 


Vn- 


+ 


bntn 


=  D. 


N 


Gauss'  elimination  method,  discussed  in  reference  4,  is  applied  to  solve 
the  system  of  equations.  This  method  affords  a  fast  and  accurate  solution  for 
matrices  containing  a  dominant  diagonal.  The  solution  of  this  matrix  gives  the 
temperature  of  each  node  in  the  system  for  the  next  future  time  step.  The 
entire  process  is  repeated  for  each  time  step  throughout  the  run,  giving  a  time 
history  of  the  temperature  at  each  node. 


Using  this  method  of  solution,  residual  errors  in  the  temperature  compu¬ 
tations  at  the  beginning  of  the  time  step  are  distributed  throughout  the  entire 
system  of  nodal  equations  and  tend  to  cancel  out  rapidly.  The  principal  advan¬ 
tage  in  using  the  implicit  method  is  a  set  of  equations  that  are  mathematically 
stable  in  time  and  distance.  Therefore,  the  magnitude  of  the  time  step  is  not 
limited  by  a  convergence  criterion.  However,  care  must  be  taken  in  selecting 
the  magnitude  of  the  time  step  in  order  to  minimize  truncation  errors  when  the 
second  derivative  of  temperature  with  respect  to  time  is  large.  A  similar 
approach  is  used  to  minimize  truncation  errors  in  distance  by  choosing  small 
node  dimensions  in  locations  where  large  second  derivatives  of  temperature  with 
respect  to  distance  are  expected. 

In  the  case  of  a  char-forming  ablative  heat  shield  where  approximately 
80  percent  of  the  heat  is  reradiated,  instability  can  arise  in  taking  large 
time  steps.  The  temperature  of  the  surface  node  can  start  oscillating  on  suc¬ 
cessive  time  steps  when  a  balance  between  the  radiation  source  and  the  heat  sink 
has  been  achieved.  Therefore,  in  ablation  problems  in  which  the  surface  node 
loses  a  large  percentage  of  heat  by  radiation,  oscillations  of  the  node  can  be 
damped  out  by  taking  small  time  steps  during-  conditions  of  high  heat  flux  and 
near  radiation  equilibrium  temperatures. 


ANALYSIS 


Figure  2  is  a  schematic  of  the  thermal  protection  system  to  be  analyzed. 

A  receding  surface  has  been  assumed  with  the  formation  of  a  residual  char  layer 
and  reaction  zone.  The  thermal  protection  system  is  composed  of  one  charring 
material  and  a  maximum  of  12  different  backup  materials  with  or  without  air  gaps. 
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The  analysis  is  such  that  the  entire  system  may  be  composed  of  noncharring 
materials.  The  thermal  properties  of  all  materials  are  temperature  dependent? 
also,  the  charring  material  properties  are  state  dependent  (fully  or  partially 
charred) . 

The  response  of  charring  ablation  heat  shields  to  a  hyperthermal  environ¬ 
ment  is  extremely  complex,  and  simplifying  assumptions  and  approximations  are 
necessary  to  afford  a  numerical  solution.  The  following  assumptions  and  approx¬ 
imations  are  utilized  in  the  equations  developed  in  this  report: 

(1)  The  material  decomposes  from  the  virgin  state  to  a  porous  char  layer 
in  the  reaction  zone. 

(2)  The  reaction  zone  can  be  defined  by  an  upper  and  lower  temperature 
limit. 

(3)  The  gas  generated  within  the  reaction  zone  is  assumed  to  pass  out  of 
the  structure  with  no  pressure  loss.  No  gas  accumulation  within  a  node  Is 
allowed. 

(4)  Local  thermal  equilibrium  is  maintained  between  the  gas  and  porous 
char  matrix. 

(5)  The  gas  undergoes  no  further  chemical  reaction  within  the  residual 
material  after  having  been  formed. 


Derivation  of  Equations 

The  equations  are  derived  for  a  moving  boundary  coordinate  system,  where 
the  front  face  is  the  moving  surface  (ref.  5)*  With  this  system,  the  ablating 
material  is  divided  into  a  fixed  number  of  nodes  of  thickness  AX  which  depends 
on  the  instantaneous  location  of  the  front  face.  The  surface  recession  is 
handled  in  a  continuous  manner,  eliminating  the  need  of  throwing  away  or  lumping 
off  of  nodes. 

The  physical  model  for  the  front  surface,  including  all  heating  terms,  is 
shown  as  follows: 
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Q 


in 


m  c  T. 
glPl  1 


p_  c  S  T 
1  P-l  1  1 


The  energy  equation  at  the  front  char  surface  is 


_d_ 

d9 


—  AX  p  c  T 
2  Tp,  1 

i  1  i 


-  5  “  hV 


1 

2 


Pnc  T 
1  Pq  1 


a  (ax) 

ae 


=  Q  +  m  c  T  +  p  S0c  T  -  m  c  T 
an  g2  p2  2  2  2  p2  2  s±  v1  1 


•  '•  /AT 

"  P1CP1S1T1  "  kl-2\£X 


(12) 


where 

S.n  ~  ^c,  blow  +  Srad  +  ^comb  ”  FecT  (T1  T°°  ) 

and 

d(AX)  d  (VL  -  S\  _  S 

d@  d§  \HP  -  l)  "  HP  -  1 

where  S  is  the  linear  surface  recession  rate  and  HP  is  the  total  number 
of  nodes  in  the  ablation  material  of  thickness  VL. 
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Rewriting  equation  (12)  in  implicit  finite  difference  form 


Q.  +  m  c  T '  -  m  c  T' 

in  S2  P2  2  gl  P1  1 


(Ti  -  Ta) 

SPlCpnTl  ~  AK  AX 

2kx  2kg 


T  *  -  T, 


p  c  s(W  ■  i-5U.  =  p  c  gV 
c  p„  -  l.oy  2  Tp  2  \ 


•  k.:'  3 


2  l  p1  1  VHP  -  1 


Then,  rearranging  and  collecting  terms  yield 


(12a) 


(121) 


The  physical  model  for  interior  points  in  the  mature  char  zone,  including 
all  heating  terms,  is  shown  in  the  following  sketch: 


•AX- 


The  energy  equation  for  interior  points  in  the  char  matrix  is 


dT.  /  • 

=  AXp.c  -  P.c  t!  ( - - — 

i  p^^  d°  ri  p±  i  yNP  -  1 


=  A  c  T '  +  k 

gi+l  pi+l  1  1 


i-l,i  \£Xj 


NP  -  1  “  f 

+  Pi+lCpi+1S  V  HP  -  1  /  Ti+1 


,  /  at\  .  , 

-  k  —  1  -me  T  -  p.c  S 

i,i+l  XJ  gi  Pi  i  i  p± 


WP  -  i  + 

g  y - £  I  ijii  (13) 

v  NP  -  1  /  i 
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Putting  equation  (lj)  in  an  implicit  finite  difference  form  yields 


(13a) 

In  the  mature  char  zone,  no  internal  gaseous  ablation  products  are  assumed 
to  form.  The  reaction  zone  is  the  source  for  the  formation  of  the  internal 
gaseous  products.  Therefore,  in  equations  (12)  and  (l3) ,  m  =  m 

gi  gi+l 

The  physical  model  for  nodes  in  the  reaction  zone  is  identical  to  sche¬ 
matic  shown  for  the  interior  nodes  in  the  char  except  for  considering  the  energy 
absorbed  in  formation  of  the  gaseous  ablation  products.  The  heat  balance 
equation  for  a  node  in  the  reaction  zone  is 


13 


Rearranging, 


The  physical  model  for  the  interface  "between  the  reaction  zone  and  virgin 
material  is  illustrated  as  follows: 
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The  heat  balance  equation  for  this  node  is 
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The  physical  model  for  an  interior  node  in  the  virgin  material  is 
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Rearranging, 


The  physical  model  for  the  last  node  in  the  ablation  material  and  first  node 
in  the  backup  structure  is 


I  f 

For  this  interface,  T.  .  =  T .  .  _ 

1,3  1,  j+1. 
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The  heat  balance  equation  for  this  node  is 


(AX .c  p.  .  +  AX  c  p, 

J  Pi^J  i.J  J+l  P±ij+1  i^J+1) 

2A0  / 


(17a) 

The  backup  structure  may  contain  up  to  a  maximum  of  12  different  materials 
with  or  without  air  gaps  between  materials.  .Therefore,  conduction  or  radiation 
and/or  convection  between  materials  is  allowed.  The  heat  balance  equations  for 
the  various  modes  of  heat  transfer  in  the  backup  structure  are  presented  in 
the  following  equations: 
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(l)  Interior  node  material: 

(t«  .  .  -  T,1  .]  IT'  .  -  T'  .] 

I  l-l, ,]  i,  ,1 J  \  i,  j  i+l,Ji  _ 


AX. 


AT 


2k.  .  .  2k.  . 

i-l,J  i,  J 


AX.  AX 

_ J_  + _  _l_ 


2k. 
i,  J 


AX. 

_ J.  I  rp  I  _  rp 

i,  3'p.,  ,  A0  |  i,  j  1,3. 
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=  P,  t. 


2k 


i+lj  j 


(18) 


Rearranging,  equation  (l8)  "becomes 
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±9  d 


T.  .  (l8a) 


(2)  First  and  last  nodes  of  two  interior  materials  with  no  gap: 


(Ti-l,,j  ~  Ti,,i)  (Ti,,i+1  ~  Ti+l,,i+l) 

AX.  AX.  "  AX.  ,  AX .  _ 

J+l 


AX,  AX,. 


.1+1 


2k.  _  .  2k.  .  2k.  2k.  _  ... 

i-l, J  i,l  i,J+l  1+1, 0+1 


P.  .c 

i,J  P 


AX.  +  P.  .  _ c  AX.  . 

i,.1  J  ^J+lpi..i+l  J+1 


2  A0 


T.*  . 

i,J 


-  T. 


i,3. 


(19) 


For  this  case,  T.'  .  =  t!  .  _ 
1,  J  1,  J+l 
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Rearranging,  equation  (19)  becomes 
1 


AX. 
_ sL 


**1-1,3  2kl,3 


T  * 

AX.  \1i-l,j 
+  - 


AX. 

_ 


AX.  AX ._  AX,.. 

_ j_  3 il  _ a±L 


+ 


.  21C.  .  2k.  •  in  1  -I  n4-l 

1-1, 1  1,1  1,1+1  1+1, J+l 


/ p.  .e  AX.  +  p.  .^nc  AX.  A 

f  1,3  pi,,i  J  1,J  1  Pi,.i+1  3  zj 

'  2  A0  ' 


T!  . 

i,l 


a: 


m.i+i  ,  M.i+i  ^i+1'J+1 


2k.  2k 

i,g+l  l+l, 1+1 


p.  .c  AX.  +  c  p.  .  .AX  ' 

^  pi„j  J  pi..i+l  1'J+1  J+1 

2  A0  i 


T.  . 
1,1 


(19a) 


(3)  First  node  of  interior  material  with  an  air  gap  between  materials: 

4  4 

T.'  .  -  T.'  ... 

1-1,J  1,1+ly 


hj(Ti-i,j  “  Ti,j+iJ  +  I  1  .  1 


ej  ei+i 


-  1 


m  *  m1 

j,.i+i  "  i+1,.1+1 

SX^,  AX... 

_ J±l_  +  _ J+l  „ 

2k.  ...  2k..  ... 

,  1,1+1  l+l, 9+1, 


p.  .  _ c  AX... 
^i.,1+1  J-1-1 

2  A9 


T'  -  T 

1  i,J+l  i,J+l 


Equation  (20)  may  be  linearized  by  using  the  approximation 

k  x  k 

T 1  =  4T  x 1  -  3T 

as  discussed  in  the  Program  Description  section. 


(20) 
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Therefore,  rearranging  and  linearizing,  equation  (20 )  becomes 


(4)  Last  node  of  an  interior  material  with  an  air  gap  between  materials: 


»■ «  I  • 

i-l,  J 


I'll  .  -  T! 

1,0  1,0+1 J 
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Rearranging  and  linearizing,  equation  (21)  becomes 


(5)  Final  node  in  backup  structure: 
(a)  Adiabatic  surface  — 


(21a) 


t'  .  -  t'. 

_ LlsL 


AX. 

. . J 


AX 


2k 


2k. 


i 


P.  .c  AX. 
P1  M  3 

- ±zj - 

2  A0 


T’.  . 


-  T. 


?-,3 


•  *1  •  £_JLV  •  . 


(22) 


Rearranging,  equation  (22)  becomes 

1 


AX.  AX. 

— J .j.  J 


2k  '  2k 

i-1 ,3  i  >3, 


AX.  AX. 

_ i_  +  _ i 


2k .  .  ’  2k.  . 

i-l, 3  1 ,3 


P.  .c  AX  A 
P1  ,  3 

— iiJ _  I  T' 

d  A6  j  i, j 


P.  .c  AX. 
P±  ,  J 
1?  J 


2  A0 


T.  . 


(22a) 
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(b)  Radiation  and/or  convection  loss  to  cabin  environment  - 


rp  *  iji  l 

i=hA  I  i>j 

AX. 

_ J_ 


AX. 


-  h  Tl  .  -  T 


env  v  i,J  env 


2k.  .  .  2k.  . 

i-l/O  i)J, 


,  4 

-  F  a  IT.  .  -  T 
env  \  i,0  envy 


p.  .c  AX. 

i,0  P,  M  0 

x)  J _  [T'  -  T 

2  A9  V  1,0  i,  0 


(23) 


Rearranging,  equation  (23)  becomes 


AX. 


\T! 


AX  .  Vi-1,  J 


h  + 


env  AX . 


—r^ —  +  F  obT?  . 
AX.  env  i,j 


2k.  .  .  2k.  . 

i-l,0  i)J, 


2k.  _  .  2k.  . 

1-1,0  i>0 


P.  .c  AX. 

1)0  P,  i  o 

-ti-U- _  I  ipj 


P.  .c  AX. 
l  -  - 


2  A^ 


i)0 


,j  Pi  i  J 

- iui -  T1 

2  A6  J  J 


1)0 


/  4  4 

-  h  T  -F  cr  3T7  .  +  T  (23a) 

env  env  env  \  i,j  env/  v  y 


Discussion  of  Assumptions 

A  brief  discussion  of  several  assumptions  and  approximations  made  in 
deriving  the  heat  balance  equations  is  now  presented. 

As  shown  in  the  Derivation  of  Equations  section,  transient  heat  con¬ 
duction,  thermal  degradation,  and  the  flow  of  the  gaseous  products  from  the 
reaction  zone  are  the  internal  thermal  transport  phenomena  of  interest. 

Several  methods  are  available  in  the  treatment  of  the  thermal  decomposition 
process,  and  they  differ  primarily  in  whether  the  chemical  decomposition 
occurs  in  a  single  plane  at  a  fixed  temperature  or  whether  a  spacially  contin¬ 
uous  decomposition  in  depth  is  assumed.  This  analysis  assumes  that  the  decom¬ 
position  from  the  virgin  to  the  char  state  occurs  in  a  reaction  zone  that  is 
defined  by  known  temperature  limits.  These  temperature  limits  are  determined 
from  thermogravimetric  test  data  for  the  particular  material  being  investigated. 
Figure  3  is  a  thermogravimetric  curve  for  typical  charring  ablation  material. 
From  this  curve,  the  rate  of  pyrolysis  m  is  calculated  by  knowing  the 
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temperature  change  of  a  particular  node  with  time,  that  is, 


Pi  "  A1 


P.'  -  P. 
i  1 


KP 

m  P  AX 

g.  ^  1  1 

1  i 

This  method  of  computing  the  gas -generation  rates  and  local  instantaneous 
density  may  be  subject  to  error  since  the  thermogravimetric  curve  of  a 
material  is  influenced  by  temperature  rise  rate  (deg/sec),  and  the  reaction 
zone  may  shift  up  and  down  the  temperature  scale.  This  error  can  be  elim¬ 
inated  by  the  use  of  an  Arrhenius  expression  of  the  form 


IS -^(p 


The  method  now  being  used  in  STAB  II  (equation  (25))  to  calculate  the  pyrolysis 
rate  is  being  investigated  to  determine  its  validity.  The  final  formulation 
of  the  pyrolysis  rate  law  must  rest  heavily  on  the  experimental  rate  data  for 
the  material  under  investigation.  The  use  of  simple  expressions  such  as 
equations  (24)  and  (25)  may  be  entirely  adequate,  depending  upon  activation 
energy  for  the  decomposition  process  and  order  of  reaction. 

The  aerodynamic  heating  input  in  the  analysis  consists  of  convective 
and  radiative  components  treated  separately.  This  distinction  is  necessary 
since  the  convective  heating  can  be  significantly  reduced  as  a  result  of  the 
injection  of  the  ablation  gases  into  the  boundary  layer,  with  generally  no 
effect  on  radiant  heating.  Reduction  in  the  convective  heating  rate  can  be 
approximated  by  the  following  expression  (ref.  6): 


q.  .  ,  =  T|m  Hm  -  H 

block  T  Wj 

Therefore, 


q  , .  =  q 

c,blow  cw 


'■block 
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However,  equation  (28)  is  unsatisfactory  for  high  blowing  rates,  since  q  block 
can  become  greater  than  4  •  An  experimental  curve  of  blocking  effectiveness 

I  ^c?blow  .  ^gHT 

yf  =  •  ac  q  fnim+.i  rm  nf  t.hp  itisqs  f.rpnqf pt*  ■npypmpt.py  — ■  y  nnn  hp  pm  — 


C,t)l0W 


as  a  function  of  the  mass  transfer  parameter 


can  be  em¬ 


ployed  to  determine  the  heating  reduction  at  high  blowing  rates.  Both  methods 
have  been  employed  in  the  STAB  II  analysis.  Equation  (28)  is  presently  in  use. 
However,  no  satisfactory  method  for  accurately  predicting  the  convective  heat 
blockage  has  been  determined. 

Another  source  of  heating  is  the  combustion  of  the  ablation  products  in  the 
boundary  layer.  Reference  7  presents  an  analysis  of  the  oxidation  of  a  carbon 
surface  and  the  resulting  combustive  heating.  The  heating  due  to  combustion  as 
derived  in  reference  7  is 

q  =  m  AH  (29) 

comb  c  c  ' 

where  AH  is  the  heat  of  combustion  per  unit  weight  of  char, 
c 

The  thermal  properties  of  the  ablation  material  are  both  temperature  and 
state  dependent  (fully  or  partially  charred) .  Figure  4  is  an  illustration  of 
the  variation  of  these  properties  with  temperature  and  state.  The  thermal  prop¬ 
erties  are  assumed  to  vary  as  follows: 


(l)  Char  zone  s  T'char) 


kc  =  f  (temp) 

c  =  f  (temp) 
P 

c 

p  =  constant 
c 


(2)  Reaction  zone  ^  T^  <  T^^  j 


p  =  f  (temp)  =  Pv  +  (Py  -  Pc) 


T.  -  T 


'abl  char 
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(3)  Virgin  zone  <  T^ 

p  =  constant 
v 

=  f  (temp) 
c  =  f  (temp) 

p 

V 

The  calculation  of  char  removal,  due  to  chemical,  thermal,  or  mechanical 
mechanism  or  a  combination  of  these  mechanisms,  has  been  examined  by  a  multitude 
of  investigators  and  numerous  correlations  exist,  depending  on  the  specific 
material  involved. 

To  provide  a  maximum  degree  of  flexibility  for  analyzing  both  ground  and 
flight  test  data  and  synthesizing  trajectories,  the  following  provisions  for 
char  removal  (surface  movement)  are  provided: 

(1)  Removal  of  char  as  a  function  of  surface  temperature. 

(2)  Removal  of  char  at  a  rate  which  is  a  function  of  time. 

As  the  char  is  removed,  the  surface  moves  with  respect  to  a  coordinate  fixed 
in  the  material.  The  distance  between  the  initial  surface  location  and  the 
char  surface  is 


ANALYSIS  VERIFICATION 


As  discussed  in  the  previous  sections,  approximations  and  assumptions 
were  made  in  the  analytical  model  to  afford  a  quick  and  accurate  solution 
in  predicting  the  thermal  response  of  a  charring  heat  shield.  These 
simplifying  assumptions  and  approximations  are  expected  to  introduce  only 
minor  errors;  however,  the  validity  of  the  analyses  and  resultant  accuracy 
can  be  judged  only  by  a  comparison  with  exact  theoretical  solutions  and 
experimental  data.  Three  examples  have  been  selected,  and  a  comparison 
of  the  STAB  II  results  with  the  theoretical  and  test  data  is  discussed  in  the 
following  paragraphs. 

An  elementary  transient  heating  example  was  chosen  to  demonstrate  the 
accuracy  and  numerical  stability  of  the  STAB  II  program.  A  steel  slab 
6  inchgs  thick  was  selected  and  assumed  to  be  at  uniform  initial  temperature 
of  46<J  R  (0°  F).  The  thermal  properties  were  considered  constant.  The 
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front  surface  was  subjected  to  a  heating  rate  of  72  Btu/ft  -sec,  and  an 
adiabatic  back  surface  was  assumed.  Figure  5  shows  a  comparison  of  the 
STAB  II  calculated  in-depth  temperatures  as  a  function  of  time  with  the 
exact  solution  taken  from  reference  8. 

To  demonstrate  the  STAB  II  solution  with  a  moving  boundary,  a  slab  with 
constant  properties,  uniform  initial  temperature,  front  surface  moving  with 
a  constant  velocity,  and  constant  surface  temperature  was  chosen.  The  exact 
solution  for  a  semi-infinite  slab  with  these  boundary  and  initial  conditions 
is  presented  in  reference  9*  Figure  6  presents  a  comparison  of  the  STAB  II 
temperature  response  with  the  exact  solutions.  As  can  be  seen  from  this 
figure,  the  two  solutions  are  not  in  agreement  for  approximately  the  first 
50  to  60  seconds  of  the  transient.  This  disagreement  is  the  result  of  the 
quasi-steady  state  assumption  made  in  the  exact  solution  analysis. 


A  calculation  was  made  to  estimate  the  induction  time  (time  at  which  ^  =  o 

is  a  good  assumption)  and  found  to  be  approximately  60  seconds,  which  is  in 
agreement  with  the  STAB  II  results. 

Finally,  to  verify  the  fully  charring  ablation  model,  an  example  of  a 
typical  charring  material  was  chosen.  (See  the  sample  problem  in  appendix  A.) 
The  charring  ablation  material  is  initially  1.6  inches  thick  with  an  adiabatic 

2 

back  surface  and  a  constant  heat  flux  of  95  Btu/ft  -sec  applied  to  the  front 
surface.  The  surface  is  assumed  to  recede  at  a  constant  velocity  of 

3*05  x  10”^in./sec.  Figure  7  presents  a  comparison  of  the  in-depth  tempera¬ 
tures  with  actual  test  results  obtained  in  an  arc  tunnel.  The  results  are  in 
good  agreement,  with  the  largest  deviations  between  calculated  and  measured 
values  occurring  for  the  thermocouple  located  at  a  depth  of  1.0  inch.  The 
disagreement  could  be  attributed  to  several  possible  errors:  thermal  property 
values,  incorrect  location  of  thermocouples,  et  cetera.  The  effect  of  varying 
the  thermal  properties  (thermal  conductivity,  specific  heat,  et  cetera)  is 
presently  being  investigated. 

Tables  I  and  II  present  the  input  and  output  data  used  for  this  example. 
Figures  8,  9  5  and  10  are  the  resulting  plot  routine  output. 

The  comparisons  between  the  computer  results  and  the  exact  solutions  and 
test  results  are  considered  satisfactory. 
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CONCLUDING  REMARKS 


An  analysis  and  a  computer  program  for  predicting  the  transient  thermal 
response  of  a  charring  ablation  thermal  protection  system  has  been  described. 
The  numerical  formulation  of  the  equations  is  such  that  an  implicit  solution 
is  obtained.  This  method  of  solution  affords  both  a  rapid  and  accurate  solu¬ 
tion  for  both  ablating  and  nonablating  type  problems . 

Provision  is  made  in  the  program  for  a  number  of  surface  boundary  condi¬ 
tions.  These  provisions  allow  efficient  use  of  the  program  for  analyzing  both 
ground  and  flight  test  data  and  trajectory  synthesis. 

The  computer  program  has  been  checked  out  with  both  exact  solutions  and 
actual  ablation  test  data.  The  numerical  results  are  in  good  agreement  with 
the  exact  solutions  and  test  data.  However,  the  analysis  depends  upon  using 
good  property  values,  and  some  effort  must  be  expended  in  obtaining  the  best 
possible  thermal  properties. 


Manned  Spacecraft  Center 

National  Aeronautics  and  Space  Administration 
Houston,  Texas,  November  1,  1965 
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APPENDIX  A 


SAMPLE  PROBLEM 


The  following  sample  problem  is  shown  to  .indicate  the  form  of  the  data 
input  and  the  program  output.  A  typical  charring  material  subjected  to  a 
constant  heating  rate  as  experienced  in  an  arc  tunnel  is  presented.  The  fol¬ 
lowing  is  a  sketch  of  the  model: 


cu  =  95  Btu/ft2-sec 

CW  ' 

T 

1.5  in. 

J_ 

0-1  in. 


The  various  material  properties  and  dimensions  are  shown  in  the  program 
output  of  Table- II.  The  insulation  is  assumed  to  be  ablation  material  for  this 
problem.  The  problem  coding  sheet  and  subsequent  data  card  listing  are  shown 
in  Table  I.  The  initial  temperature  of  the  structure  was  assumed  uniform  and 
equal  to  530°  R  (70°  F).  Figures  8,  9  5  and  10  are  the  output  data  obtained 
from  the  plot  routines. 
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APPENDIX  B 


PROGRAM  USAGE  INSTRUCTIONS 


IBM  'JO^k-fk-O  program  F021,  standard  ablation  program,  designated  STAB  II, 
is  designed  to  evaluate  the  transient  thermal  performance  of  a  charring  ablation 
heat  protection  system.  The  program  considers  ojie  ablating  material  and  up  to 
12  different  materials  in  the  supporting  backup  structure.  A  maximum  of 
50  nodes  may  be  considered  in  the  ablation  material,  and  a  maximum  of  10  nodes 
per  material  is  allowed  for  each  backup  structure  material.  Air  gaps  can  be 
considered  between  successive  materials  in  the  backup,  thus  allowing  for 
both  radiative  and/or  convective  heat  transfer  between  materials.  The  heat  loss 
to  the  cabin  environment  from  the  backup  structure  can  be  accomplished  by  both 
radiation  and/or  convection,  or  an  adiabatic  backface  surface  may  be  prescribed. 

Unless  otherwise  specified,  the  input  problem  data  are  in  "floating 
point"  form  (E12.8  format)  and  must  end  in  columns  12,  24,  36,  48,  60,  and  72. 

It  is  suggested  that  each  floating  point  number  have  a  sign,  a  two-digit 
exponent,  and  a  decimal  point.  For  example,  the  number  145- 23  can  be  written 
as  +1.4523  +02,  +145.23  +00,  or  +.14523  +03. 


Input  Nomenclature 

The  nomenclature  used  in  the  problem  data  input  is  as  follows: 

NCASE  number  of  problems  to  be  run  successively 

HEAD  any  72  alphabetical  and/or  numerical  characters 

TITLE  control  card  for  reading  in  new  input  for  successive 
problems 

1.  blank  card  —  new  data  will  be  read  in 

2.  Six  asterisks  in  columns  1  to  6.  Skip  to 

next  read  statement 

TLIM  time  limit  of  problem,  sec 
TINT  starting  time  of  problem,  sec 

NPTT  number  of  points  in  time— step  table  (the  minimum 

value  of  NPTT  is  2) 

NPL/T  output  plot  control 

=1  plot  routine  will  be  used 

=0  plot  routine  will  be  ignored 
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TTABLE 


DELTT 

IPRC 

FC0KV 

FRAD 

TABL 

TCHAR 

TREC 

RH0V 

RH0C 

FBL0W 

EMV 

EMC 

H300 

VL 

HV 

VPT 


=0 

=1 

FV 

TV 


time  in  time— step  table,  sec 

time  step  to  be  used  for  each  calculation  —  starting 
at  time  TTABLE,  sec 

variable  print  frequency  in  TTABEE  table;  that  is,  if 
DELTT  =  1. 0  and  IPRC  =  10,  the  output  will  be  printed 
at  10— second  intervals 

factor  to  correct  convective  heating  rate  for  various 
body  locations 

factor  to  correct  radiative  heating  for  various  body 
locations 

temperature  at  which  ablation  starts,  °R 

temperature  at  which  ablation  stops,  °R 

surface  temperature,  °R,  or  time  at  which  char  removal 
is  to  start,  sec 

density  of  virgin  ablation  material,  lb/ft^ 
density  of  mature  char  material,  lb/ft ^ 

blowing  efficiency  of  ablation  gases  in  reducing  convective 
heating 

emissivity  of  virgin  ablation  material 

emissivity  of  charred  ablation  material 

enthalpy  of  air  at  300°  K,  129.  06  Btu/lbm 

initial  thickness  of  virgin  ablation  material,  in. 

heat  of  degradation  of  virgin  material,  Btu/lbm 

test  to  determine  if  the  reaction  zone  and  char  zone 
thermal  properties  are  irreversible  with  temperature 

properties  are  irreversible  and  equal  to  the  value  at 
the  maximum  individual  node  temperature  (this  is  the 
recommended  value  for  VPT) 

properties  are  reversible 

view  factor  for  external  environment 

sink  temperature  of  external  environment,  °R 
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CHARE  thermal  conductivity  of  material  at  TCHAR,  Btu/ft— hr-°R 

CHARC  specific  heat  of  material  at  TCHAR,  Btu/lhm-°R 

ABLK  thermal  conductivity  of  material  at  TABL,  Btu/ft— hr— °R 

ABLC  specific  heat  of  material  at  TABL,  Btu/lhm-°R 

RP  number  of  node  points  in  ablation  material 

RKC  number  of  points  in  char  thermal  conductivity  —  temper¬ 

ature  table 

RCPC  number  of  points  in  char  specific  heat  —  temperature 

table 

REV  number  of  points  in  virgin  thermal  conductivity  - 

temperature  table 

RCPV  number  of  points  in  virgin  specific  heat  —  temperature 

table 

HREC  number  of  points  in  surface  recession  —  temperature  or 

time  table 

TEC  temperature  values  in  char  thermal  conductivity  - 

temperature  table,  °R 

XKC  thermal  conductivity  values  in  char  thermal  conductivity 

temperature  table,  Btu/ft— hr— °R 

TCPC  temperature  values  in  char  specific  heat  —  temperature 

table,  °R 

CPC  specific  heat  values  in  char  specific  heat  —  temperature 

table,  Btu/lbm-cR 

TEV  temperature  values  in  virgin  thermal  conductivity  — 

temperature  table,  °R 

XKV  thermal  conductivity  values  in  virgin  thermal  conduc¬ 

tivity  —  temperature  table,  Btu/ft— hr— °R 

TCPV  temperature  values  in  virgin  specific  heat  —  temperature 

table,  °R 

CPV  specific  heat  values  in  virgin  specific  heat  temperature 

table,  Btu/lbm— °R 

TS  temperature, °R,  or  time,  sec,  values  in  the  surface 

recession  table 


SR 

KTRAPT 

TIME 

QC0N 

QRAD 

VEL 

KMB 

NEBS 

BL 

XNPM 

NKPB 

KCPB 

XIDNT 

TXK 

XK 

TCP 

CPX 

RH0BX 

XBM 

E1MFB 


surface  recession  values  in  the  surface  recession  — 
temperature  or  time  table,  in. /sec 

number  of  time  points  in  the  trajectory  input  table 

the  array  of  (WTRAPT)  trajectory  time  values,  sec 

the  corresponding  array  of  cold  wall  convective  heating 
2 

rates,  Btu/ft  -sec 

the  corresponding  array  of  radiative  heating  rates, 
Btu/ft*'-  sec 

the  corresponding  array  of  flight  velocity,  ft/sec 

number  of  materials  in  backup  structure 

total  number  of  node  points  in  backup  structure 

total  thickness  of  backup  structure,  in. 

number  of  nodes  in  each  individual  material  in  backup 
structure 

number  of  points  in  each  individual  backup  structure 
material  thermal  conductivity  -  temperature  table 

number  of  points  in  each  individual  backup  structure 
material  specific  heat  —  temperature  table 

any  72  alphanumeric  characters  used  to  describe  each 
individual  material  in  the  backup  structure 

temperature  values  in  backup  material  thermal  conductivity 
temperature  table,  °R 

thermal  conductivity  values  in  backup  material  thermal 
conductivity  —  temperature  table,  Btu/ft— hr— °R 

temperature  values  in  backup  material  specific  heat  — 
temperature  tables,  °R 

specific  heat  values  in  backup  material  specific  heat  — 
temperature  tables,  Btu/lbm-°R 

density  of  individual  materials  in  backup,  Ib/fu 

thickness  of  individual  materials  in  backup,  in. 

emissivity  of  front  surface  of  each  material  in  backup 
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EMBB  emissivity  of  back  surface  of  each  material  in  backup 

H  film  coefficient  between  adjacent  materials  in  backup, 

Btu/ft^-hr-°R 

GAPX  width  of  gap  between  adjacent  materials  in  backup,  in. 

FTEST,  tests  to  determine  the  mode  of  heat  transfer  between 
BTEST  materials  for  the  front  and  backface  of  each  material 

respectively 

=0  conduction  only  between  materials 

=+l  convective  heat  transfer  only 

=-l  radiation  only  or  radiation  and  convection  heat  transfer 

TEW  temperature  of  interior  cabin  environment,  °R 

HEW  film  coefficient  to  interior  cabin  environment, 

Btu/ft2-hr-°R 

FENV  view  factor  and  emissivity  product  for  radiative  heat 

transfer  to  cabin  interior 


QL0SS 


boundary  condition  between  last  node  of  the  backup 
structure  and  cabin  environment 


=0  adiabatic  surfaces 

=+l  radiation  and/or  convective  loss 

TEST2  determines  the  proper  heat  shield  initial  temperature 
distribution 

=0  constant,  uniform  initial  temperature  distribution 

=-l  arbitrary  initial  temperature  distribution 

=+l  linear  temperature  distribution 

TEMPI  temperature  to  be  used  when  constant  temperature 
distribution  option  is  used,  R 

TX0  initial  temperature  at  front  surface  of  heat  shield  to 

be  used  in  computing  initial  linear  temperature 
gradient ,  R 

TEMDI  arbitrary  temperature  distribution  values,  to  be  used  only 
if  TEST2  is  negative,  °R 


NHP  number  of  points  in  enthalpy  —  temperature  curve  fit 

HX  enthalpy  values  in  enthalpy  —  temperature  table, 

Btu/lbm 

TW  corresponding  temperature  values  in  enthalpy  —  temperature 

table,  R 


Input  Data  Card  Preparation 

The  input  data  are  given  in  the  following  order.  Each  number  in  the  fol¬ 
lowing  listing  refers  to  a  separate  record  and  must  begin  on  a  new  data  card. 
The  input  data  have  been  grouped,  where  possible,  into  various  sections  dealing 
with  a  particular  part  of  the  input,  that  is,  ablation  material  properties, 
trajectory  data,  backup  structure,  et  cetera.  This  grouping  permits  the  use  of 
a  minimum  number  of  input  cards  for  running  successive  problems.  The  title 
card  as  described  in  the  input  nomenclature  controls  the  input  for  successive 
problems . 


1.  The  first  data  card  contains  the  value  of  NCASE.  NCASE  is  an  integer 
(15  format)  and  must  end  in  column  5<  This  card  tells  how  many  problems  are  to 
be  run  and  is  entered  only  once  at  the  start  of  the  data  deck. 

2.  Columns  1  to  72  of  the  second  data  card  contain  any  title  or  identi¬ 
fication  information  desired;  any  alphanumeric  character  may  be  used.  This 
card  is  printed  at  the  top  of  the  first  page  of  the  output.  This  card  must  be 
included  in  all  successive  problems  to  be  run. 

(a)  Problem  time  section 

3.  TITLE  card  —  if  blank,  cards  4  and  5  must  be  submitted.  If  six  aster¬ 
isks  are  punched  in  columns  1  to  6,  skip  to  record  number  6. 

4.  This  record  contains,  in  the  following  order,  TLIM,  TINT,  NPTT,  and 
NPL0T.  TLIM  and  TINT  are  entered  as  floating-point  numbers  and  must  end  in 
columns  12  and  24.  NPTT  and  NPL0T  are  integers  entered  with  an  15  format  and 
must  end  in  column  30  and  35* 

5.  Start  entering  the  values  of  TTABLE,  DELTT,  IPRC.  TTABLE  and  DELTT 
are  floating-point  numbers  and  must  end  in  columns  12  and  24.  IPRC  is  entered 
as  integer  with  an  15  format  and  must  end  in  column  30.  Use  as  many  cards  as 
required  to  enter  NPTT  values. 

(b)  Heating  rate  factors  section 

6.  TITLE  card  -  if  blank,  card  7  must  be  submitted.  If  six  asterisks  are 
punched  in  columns  1  to  6,  skip  to  record  number  8. 


38 


7.  Enter  the  FC0NV  and  FRAD.  These  numbers  are  entered  as  floating-point 
numbers  and  must  end  in  columns  12  and  24. 

(c)  Ablation  material  section 

8.  TITLE  card  _  if  blank,  cards  9  to  18  must  be  submitted.  If  six  aster¬ 
isks  are  punched  in  columns  1  to  6,  skip  to  record  number  19 . 

9.  HEADNG  card  —  any  alphanumeric  characters  in  columns  1  to  72.  Records 
9  to  18  contain  input  data  for  the  ablation  material. 

10.  Enter  TAEL,  TCHAR,  TREC,  RH0V,  RH0C,  and  FBL0W.  These  numbers  are 
entered  as  floating-point  numbers  (6E12.8  format)  and  must  end  in  columns  12, 

24,  56,  48,  60,  and  72. 

11.  Enter  EMV,  EMC,  H300,  VL,  HV,  and  VPT.  Use  the  same  format  as  card  10. 

12.  Enter  FV,  TV,  CHARK,  CHARC,  ABLE,  and  ABLC .  Use  the  same  format  as 
card  10. 


13.  This  card  contains,  in  the  following  order,  UP,  MC,  NCPC,  NKV,  NCPV, 
and  NREC.  These  numbers  are  fixed-point  integers  and  must  end  in  columns  5? 

10,  15,  20,  25,  and  30.  An  15  format  is  used  to  read  in  these  numbers. 

14.  Start  entering  the  curve  of  TKC  versus  XKC,  with  the  values  of  TKC 
ending  in  columns  12,  38,  and  60.  The  corresponding  values  of  XKC  must  end  in 
columns  24,  48,  and  72;  for  example,  three  TKC-XKC  points  are  contained  on  one 
card.  The  numbers  are  entered  as  floating-point  numbers.  Use  as  many  cards  as 
required  to  enter  NKC  points  on  the  curve. 

15*  Start  entering  the  curve  of  TCPC  versus  CPC  with  the  values  of  TCPC, 
ending  in  columns  12,  36,  and  60.  The  corresponding  values  of  CPC  must  end 
in  columns  24,  48,  and  72;  for  example,  three  TCPC-CPC  points  are  contained  on 
one  card.  The  numbers  are  entered  as  floating-point  numbers.  Use  as  many 
cards  as  required  to  enter  the  NCPC  points  on  the  curve. 

l6.  Start  entering  the  curve  of  TKV  versus  XKV  with  the  values  of  TKV 
ending  in  columns  12,  38,  and  60.  •  The  corresponding  values  of  XKV  must  end  in 
columns  24,  48,  and  72;  for  example,  three  TKV-XKV  points  are  contained  on  one 
card.  The  numbers  are  entered  as  floating  point.  Use  as  many  cards  as  re¬ 
quired  to  enter  the  NKV  points  on  the  curve. 

17-  Start  entering  the  curve  of  TCPV  versus  CPV  with  the  values  of  TCPV, 
ending  in  columns  12,  38,  and  60.  The  corresponding  values  of  CPV  must  end  in 
columns  24,  48,  and  72;  for  example,  three  TCPV-CPV  points  are  contained  on  one 
card.  The  numbers  are  entered  as  floating  point.  Use  as  many  cards  as  re¬ 
quired  to  enter  NCPV  points  on  the  curve. 

18.  Start  entering  the  curve  of  TS  versus  SR  with  the  values  of  TX,  ending 
in  columns  12,  36,  and  60.  The  corresponding  values  of  SR  must  end  in  columns 
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24,  48,  and  72;  for  example,  three  TS-SR  points  are  contained  on  one  card. 

The  numbers  are  entered  as  floating  point.  Use  as  many  cards  as  required  to 
enter  NREC  points  on  the  curve. 

(d)  Trajectory  data  section 

19.  TITLE  card  -  if  blank,  cards  20  to  22  must  be  submitted;  if  six 
asterisks  are  punched  in  columns  1  to  6,  skip  to  record  number  23. 

20.  HEADNG  card  -  any  alphanumeric  characters  in  columns  1  to  72.  Records 
21  and  22  contain  trajectory  input  data. 

21.  Enter  NTRAPT.  This  number  is  an  integer  and  must  end  in  column  5 • 

An  15  format  is  used  to  read  in  this  number. 

22.  Start  entering  the  trajectory  data  in  the  following  order:  TIME, 

QC0N,  QRAD,  VEL.  These  values  are  entered  as  floating-point  numbers  and  must 
end  in  columns  12,  24,  36,  and  48.  There  are  four  trajectory  data  points  on 
one  card.  Use  as  many  cards  as  required  to  enter  NTRAPT  points  in  the  tra¬ 
jectory. 

(e)  Backup  structure  section 

23.  TITLE  card  -  if  blank,  cards  24  to  31  must  be  submitted;  if  six 
asterisks  are  punched  in  columns  1  to  6,  skip  to  record  number  32. 

24.  Enter  NMB,  NPBS,  and  BL.  These  three  values  must  end  in  columns  5? 

10,  and  24.  NMB  and  NPBS  are  integers  and  are  read  in  under  an  15  format.  BL 
is  a  floating-point  number. 

25.  Enter  the  values  of  XNPM.  XNPM  is  in  floating-point  form  and  must 
end  in  columns  12,  24,  36,  48,  60,  and  72.  Use  as  many  cards  as  required  to 
enter  NMB  points. 

26.  Enter  the  values  of  NKPB  and  NCPB.  These  numbers  are  integers  and 
NKPB  must  end  in  columns  55  15?  25,  35 >  and  45;  and  the  corresponding  values  of 
NCPB  must  end  in  columns  10,  20,  30,  40,  and  50.  An  15  format  is  used  to  read 
these  values.  Five  NKPB-NCPB  values  are  contained  on  one  card.  Use  as  many 
cards  as  are  required  to  enter  NMB  points. 

27.  XIDNT  card  -  any  alphanumeric  characters  in  columns  1  to  72.  This 
card  contains  a  description  of  each  backup  material. 

28.  Start  entering  the  curve  of  TXK  versus  XK  with  the  values  of  TXK, 
ending  in  columns  12,  36,  and  60„  The  corresponding  values  of  XK  must  end  in 
columns  24,  48,  and  72;  for  example,  three  TXK-XK  points  are  contained  on  one 
card.  The  numbers  are  entered  as  floating  point.  Use  as  many  cards  as  re¬ 
quired  to  enter  NKPB  points  on  the  curve. 
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29.  Start  entering  the  curve  of  TCP  versus  CPX  with  the  values  of  TCP, 
ending  in  columns  12,  36,  and  60.  The  corresponding  values  of  CPX  must  end 
in  columns  24,  48,  and  72;  for  example,  three  TCP-CPX  points  are  contained 
on  one  card.  The  numbers  are  entered  as  floating  point.  Use  as  many  cards  as 
required  to  enter  NCPB  points  on  the  curve.  Repeat  records  27,  28,  and  29 
until  the  properties  for  NMB  materials  have  been  entered.  The  maximum  number 
for  RMB  is  12. 

30.  Start  entering  the  following  values  in  order:  RH0BX,  XBM,  EMFB,  and 
EMBB.  These  values  are  entered  as  floating-point  numbers  (6E12.8  format) 
and  must  end  in  columns  12,  24,  36,  48,  60,  and  72.  Use  as  many  cards  as 
required  to  enter  EMB  points. 

31.  Start  entering  the  following  values  in  order:  H,  GAPX,  FTEST,  and 
BTEST.  These  values  are  entered  as  floating-point  numbers  (6E12. 8  format) 
and  must  end  in  columns  12,  24,  36,  48,  60,  and  72.  Use  as  many  cards  as 
required  to  enter  EMB  points. 

(f)  Interior  environment  section 

32.  TITLE  card  —  if  blank,  cards  33  and  34  must  be  submitted;  if  six 
asterisks  are  punched  in  columns  1  to  6,  skip  to  record  number  35. 

33*  HEADNG  card  -  any  alphanumeric  characters  in  columns  1  to  72. 

Record  35  contains  properties  of  environment. 

34.  Enter  the  following:  TEEV,  HEEV,  FEEV,  and  QL0SS.  The  values  are 
entered  as  floating-point  numbers  and  must  end  in  columns  12,  24,  36,  and  48. 

(g)  Initial  temperature  section 

35*  TITLE  card  —  if  blank,  records  36  and  37  must  be  submitted;  if  six 
asterisks  are  punched  in  columns  1  to  6,  skip  to  record  39. 

36.  HEADEG  card  —  any  alphanumeric  characters  in  columns  1  to  72. 

Records  37  and  38  contain  initial  temperature  distribution  input. 

37.  Enter  TEST2,  TEMPI,  and  TX0.  These  values  are  entered  as  floating¬ 
point  numbers  and  must  end  in  columns  12,  24,  and  3 6. 

NOTE:  If  TEST2  is  a  negative  number,  record  38  must  be  submitted;  other¬ 
wise,  skip  to  record  39 '. 

38.  Enter  the  arbitrary  temperature  distribution  values,  TEMDI.  These 
values  are  entered  as  floating  points  with  a  6E12. 8  format.  Use  as  many 
cards  as  required  to  enter  EP  plus  NFBS  node  points. 

(h)  Enthalpy  -  temperature  section 


4l 


39 •  TITLE  card  —  if  blank,  records  40  and  4l  must  be  submitted;  if  six 
asterisks  are  punched  in  columns  1  to  6,  this  is  the  last  data  card  in  the 
problem  input. 

40.  Enter  NHP.  This  value  is  an  integer  and  must  end  in  column  5*  An 
15  format  is  used  to  read  in  this  number. 

41.  Start  entering  the  curve  of  HX  versus  TW  with  the  value  of  HX  ending 
in  columns  12,  5 6,  and  60.  The  corresponding  values  of  TW  must  end  in  columns 
24,  48,  and  72;  for  example,  three  HX-TW  points  are  contained  on  one  card. 

The  numbers  are  entered  as  floating  points.  Use  as  many  cards  as  required  to 
enter  NHP  points  on  the  curve.  Record  4l  consists  of  the  last  data  cards 
required  as  input  for  a  problem. 

As  many  successive  problems  as  desired  may  be  run  at  one  time  by  proper 
input  preparation.  STAB  II  has  been  designed  to  save  all  input  information 
until  it  is  changed  by  new  input  data.  Therefore,  the  use  of  the  TITLE 
control  card  is  very  important  when  running  more  than  one  problem  and  using 
the  input  data  of  previous  problems.  As  shown,  each  input  section  starts 
with  a  TITLE  control  card  for  determining  whether  new  input  data  are  to  be 
used.  If  any  data  are  changed  within  a  section,  then  all  data  cards  required 
for  that  section  must  be  submitted. 

STAB  II  can  also  be  used  for  solving  one-dimensional  transient  heat- 
conduction  problems  of  nonablating  materials.  The  following  input  parameters 
must  be  adhered  to: 

1.  TABL  must  be  greater  than  the  maximum  temperature  expected  during  the 
calculation.  Also,  TABL  >  TCHAR  >  TREC. 

2.  The  ablation  material  must  be  considered  to  be  the  first  material 
in  the  structure  for  calculation  purposes. 

5-  The  virgin  and  char  properties  must  be  inputed  as  described  above 
but  can  have  the  same  values;  that  is,  XfCV  =  XKC,  CPC  =  CPV,  RH0V  =  RH0C, 
et  cetera. 

The  following  dimensional  statements  and  program  limitations  should  not 
be  violated  when  preparing  the  input  described  above  for  ablating  and  non- 
ablating  structure: 

1.  All  property  tables  can  have  a  maximum  of  20  points  (i. e. ,  a  tempera¬ 
ture  and  specific  heat  value  constitute  one  point). 

2.  The  surface  recession  table  can  have  a  maximum  of  50  points  (TS  and 
SR  constitute  one  point). 

3.  The  trajectory  table  can  have  a  maximum  of  300  points  (TIME,  QC0N, 
QRAD,  and  VEL  constitute  one  point). 
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4.  The  ablation  material  can  be  broken  into  a  maximum  of  50  nodes.  The 
backup  structure  can  consist  of  up  to  12  different  materials  with  a  maximum  of 
10  nodes  per  material. 

5-  A  minimum  of  three  nodes  per  material  (ablation  or  backup)  must  be 
specified. 

6.  A  minimum  of  two  materials  must  be  specified  (ablation  material  and 
one  backup  structure  material). 

7.  Pure  conduction  only  is  allowed  between  the  ablation  material  and  the 
first  material  in  the  backup. 

8.  If  any  data  input  is  changed  in  the  Ablation  Material  Section  on 
successive  problems,  the  Ablation  Material  Section  data  cards  plus  the 
Initial  Temperature  Section  data  cards  must  be  submitted. 


Program  Output  Information 

The  computed  results  are  available  in  two  forms  of  output:  tabular  and 
plot  outputs.  The  tabular  output  presents  the  computed  results  in  block 
type  form  for  each  computation  step  as  controlled  by  the  print  count  control 
number.  As  discussed  in  the  preparation  of  input  data,  both  the  computational 
time  step  and  print  control  can  be  varied  throughout  the  running  of  a  problem. 
Therefore,  excessive  printed  output  is  avoided,  and  there  is  a  considerable 
savings  in  actual  machine  computation  time.  The  plot  outputs  are  printed 
and  plotted  only  when  the  entire  set  of  problems  to  be  run  are  completed. 

Tabular  output.  —  The  program  prints  a  listing  of  the  data  input  para¬ 
meters  for  identification  of  the  problem  and  ease  in  identifying  any  input 
mistakes.  For  stacked  problems,  the  program  prints  only  that  input  information 
that  is  changed  from  the  previous  problem.  The  following  calculated  problem 
output  is  printed: 

1.  Time,  sec 

p 

2.  Cold  wall  convective  heating  rate  without  blowing,  Btu/ft  -sec 

3.  Radiative  heating  rate,  Btu/ft^- sec 

4.  Velocity,  ft/sec 

5.  Gas  a Hation  rate,  lbm/ft^-hr 

6.  Char  ablation  rate,  lbm/ft  -hr 

7.  Total  ablation  rate,  lbm/ft^-hr 

8.  Surface  recession  depth  from  original  surface,  in. 
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9-  Hot  wall  convective  heating  rate  without  "blowing,  Btu/ft  -sec 

10.  Temperature  distribution  in  ablation  material,  °R 

11.  Temperature  distribution  in  backup  structure,  °R 

The  temperatures  printed  for  the  ablation  material  are  for  fixed  distances 
from  the  original  surface.  These  distances  are  calculated  from  the  initial 
ablation  material  thickness  and  number  of  nodes  in  ablation  material.  For  ex¬ 
ample 

let 


VL  =  1.0  in. 


HP  =  11 


then 


AX  = 


VL 

HP  -  1 


0.1 


The  temperatures  will  be  printed  for  X  distances  of  0,  0.1,  0.2,  0.3, 
et  cetera,  from  the  original  surface  until  the  surface  has  receded  beyond  these 
fixed  distances  at  which  time  the  node  no  longer  exists  and  is  dropped  from 
the  printout.  This  is  illustrated  in  the  following  way:  let  surface 
recession  =  0.26  inch.  The  first  temperature  printed  then  is  the  surface 
temperature  of  the  material,  located  0. 26  inch  from  the  original  material  sur¬ 
face.  The  following  printed  ablation  material  temperatures  are  for  X  dis¬ 
tances  of  0.3,  0.4,  0.5,  . ..,  1.0  inch. 

The  format  for  the  temperature  distribution  printout  is  El6. 5  with  six 
temperatures  printed  per  line. 

Plot  output.  —  The  plot  output  gives  the  following  ablative  material  per¬ 
formance  parameters  as  a  function  of  time: 

1.  Surface  depth,  in. 

2.  Bondline  temperature  between  ablator  and  backup  structure,  °R 

3-  Two  selected  isotherm  depths 

These  values  are  also  printed  in  tabular  form  for  ease  in  checking  and 
replotting  of  the  results.  The  plotted  curves  contain  all  maximum  and  mini¬ 
mum  values  of  the  parameters. 
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APPENDIX  C 

PROGRAM  IN  FORTRAN  STATEMENTS 


^5 


UOO<JOOO<J<J 


SIBFTC  MAIN 

STRUCTURES  AND  MECHANICS  PIVISION 
THERMO-STRUCTURES  BRANCH 
THERMAL  PROTECTION  SYSTEMS  SECTION 


THIS  PROGRAM  DETERMINES  THE  PERFORMANCE  OF  A  CHARRING  APlATOR 
ANALYSIS  AND  PROGRAM  DEVELOPED  Py  DONALD  M.  CURRY  *  ES32 


C 


DIMENSION  ESAVE1 ( 3) f ESAVF?  <3) ,ESAVE3(3> 

DIMENSION  TITLE <121 ,HFADNG(12) , XlDNT ( 12 , 1 2 ) ,TKC(20) »XKC(20) , 
lc£>C<20)  »TKV(P0>  ,XKV(pO)  .TCPV(PO)  ,CPV<20) , TIME ( 30n >. GCON ( 300 ) . 

20RAD ( 300 ) »VEL(300) » XnPM < 1 2 ) .NKPP ( 1 2 ) , NCPR ( 12 ) .TXk ( 20 , 12 ) , XK { 20 , 12 > 

3.  TCP ( 20 • 12) ,CPX<20,12) , RHOBX ( 1 2 ) , XBM < 1 2 ) ,EMFB(12) »EM8R(12) ,HXX()2) 

4 ,  GAPX ( 12 ) »FTFST (12)»fiTEST(12) .TEmDI ( 200 ) , TXl (200) »TX2(200) , 

5TX2T (10*12) >  TIIL1  ( 200 ) , TUL2 (200)*HX(50) »Tw( 50 > , IP ( SO ) , IR1 ( 50 ) , 
6lR2(50) #TUL(SO) * IFM ( SO ) , TY ( 200 ) , A ( 200 ) .P ( POO ) * C ( POO ) ,0(200) , 

7R ( 50 ) , RHo (50 ) ,CP(50) ,DXP( 12) ,XKP(1 0,12) , CPB( 10, 12) ,XMDG(50) , 
flYK (50) ,AB(10,1 2) , BB ( 10, 12) »CR(10,12) ,DB(10,12) , SR (10, 12) , 
9RP1(10,12) *RB 2(10,12) ,H(12) ,S(50) , NPM ( 12 ) 

DIMENSION  TTUL ( 50 ) *RHOY 1 ( SO ) , PH0y2 ( 50 ) .DrHO ( 50 ) ,TCPC ( 20 ) 

DIMENSION  TImFP (300),PRES(300),Xc(50) ,TX2C(50) , XV ( 50 ) , XDV t 50 ) 
DIMENSION  TS ( SO ) , SR ( 50 ) 

DIMENSION  TTARLE(20). DELTT (20) *IPPC(20) 

DIMENSION  ASAVE1 (3) , ASAVF2 < 3 ) , ASAVE3 ( 3 > ,BSAVE1 (3) ,BSAVE2(3) , 
1HSAVE3(3) ,CSAVE1(3) ,CSAVE2<3) ,CSAVE3(3) ,HEAD( 12) , 

1DSAVEK3)  ,DSAVE2 (  3 )  ,0SAVF3(3) 

DIMENSION  XRA ( 30 ) , Y A ( 30 ) 


COMMON  TkC , XkC , TCPC , CPC , TKV , XK V , TCPV , CPV, XNPM , RHoRX , XRM, EMBR , 
1EMFP, NKPb,NCPR, TXK, Xk , TCP, CPX, NPM* GAPX, FTFST . RTFsT , TEMDI , TXl , 
2TX2 , TX2T , TUL , TUL1 , TUl.2, IR , IR1 , IRP , A , B,C ,D, S , R , AR , RB, CR , DB , SR , 
3RR1.RR2,TY,RHOY1.RH0Y2.XMDG,RH0,CP*YK,XKR,CPB,DXr,DT,XL0ST. 
4TABL,TCHAR,TREC,RHOV,RHOC,FBLDW.fMV,EMC,H30O,NKC,NCPC,NKV,NCPV, 
5NP , NMB , NPPS , NPF , TFSTp , TFMPI , TXO , TFNV , HENv , FFNV , GLOSS , TLlM , TINT 
COMMON  Il*I2,T3,I4»l5,I6»0lN,TNT,DX»XMT,TL*VL» BL,DMP, FRR 1 ,ERR2, 
1FPR3,ERR4,HV,VPT,CHARK,CHARC. ABLk*ABLC*XmDC,H 


C 

3000  FORMAT ( 12A6 ) 

3001  FORMAT ( IX, 12AS) 

3002  F0RMAT(6E12.8) 

3003  FORMAT (615 ) 

3004  FORMAT (1X5) 

3005  FORMAT ( 215 ) 

3007  FORMAT (2I5»lE14.8) 

3008  FORMAT (///1X.12A6) 

3009  FORMAT (1H1»1X,12A6) 

3010  FORMAT (4E12»8) 

3011  FORMAT (2E12, 8, 16 , I5,E1 3. 8 ,E12 .8 ) 

3012  F0RMAT(PE12.8,I6> 

DATA  PRV0US/0545u545454S4/ 

RFWIND  11 

ST0P=9999. 

READ (5,3003) NCASE 
LPLOTrO 


JCNT=0 


AOOOO 

A0010 

A0020 

A0030 


AO  040 
A0050 
A0060 
A0070 


A0080 

A0090 

AOlOO 

A0110 

A0120 

A0130 

A0140 

A0150 

A0160 

A0170 

A0180 

A0190 

A0200 

A0210 

A0220 

A0230 

A0240 

A0250 


A0260 

A0270 

A0280 

A0290 

A0300 

A0310 

A0320 

A0330 

A0340 

A0350 

A0360 

A0370 

A0380 

A0390 

A0400 

A0410 

A0420 

A0430 


A  0440 
A0450 
A  0460 
A0470 
A0480 
A0490 
A0500 
A0S10 
A0520 
A0530 
A0540 
A0550 
A0560 
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50  NK  =  1  A0570 

11=2  *0580 

12=2  A0590 

13=2  A0600 

14=2  A0610 

15=2  A0620 

I 5=2  A0630 

117=2  A0640 

INT=1  A0650 

XLOST=0 • 0  A0660 

XMT=0.0  A0670 

XMDT=0.0  A06R0 

FRR1=0.0  A0690 

FRR2=0.0  A0700 

FPR3=0 • 0  A0710 

FRR4=f).0  A0720 

ICT=0  A0730 

ICONT=0  A0740 

XMDC=0.0  A0750 

NKP=1  A0760 

XL  STV=0 • 0  *0770 

NRS=2  A0780 

FRR5=0 . 0  A0790 

TPCT=0  A0800 

ICTP=0  A0810 

IPL0T=1  A0820 

NXA=1  A0830 

NXB=1  A0840 

NXC=1  A0850 

NXO=L  A0880 

8AVY3=-loO.  A0870 

9AVY4X=“100 .  A0880 

9X0=0. 0  A0890 

9DOT=0.0  A0900 

A0910 

GFNFRAL  TITLF  OF  PROBLEM  A0920 

100  READ ( 5, 3000 )  ( HEAD (K ) ,K=1 . 12 )  A0930 

WRITE (6,3009)  ( HEAD (K ) , K=1 . 12 )  A0940 

LPLOT=LPLOT+l  A0950 

WRITE  (11) (HEAD(I) ,1=1,12)  A0960 

WRITE (6,110)  A0970 

1)0  FORMAT(//lX»l 1 HINPUT  OATA.//)  A0980 

HEAD ( 5, 3000 )  (TITLE(L) ,L=1 ,12)  A0990 

IF (TITLE ( U .F0.PRV0U9)  GO  TO  150  A1000 

RFADCS,3011)TLIM,TINt,NPTT,NPLOT,DMP,TDMP  A1010 

RFAD( 8, 3012 )  (TTARLEt I ) , DFLTT ( I ) , IPRC ( I ) , 1=1 ,NPTt)  A1020 

T=TINT  A1030 

OTS=DFLTT(l)  A1040 

DT=DELTT (1 )/3ft00.0  A1050 

WRITE (6, 120 )  TLIM.TINT.NPTT  A1060 

120  FORMAT(1HO,11HTIME  LIMIT=,1PF 10. 4*4X»!3H INITIAL  TlMEr, 1PE10 .4 , 4X . 8  A1070 

lHMPTT= , I 4 )  A1080 

WRITE (6 , 122 )  A1090 

122  FORMAT(//8X,4HTIME,10X,9HTIMF  STfP,6X, 13HPRINT  CONTROL)  A1100 

WRITE(6, 124)  (TTAPLE(T) ,DFLTT(I) ,IPRC(I) ,I=1,NPTT)  AlilO 

124  FORMAT (5X,lPEl0.4,6Xf) PF 10.4,9X,I4)  A1120 

C  A1130 
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c  location  factors  for  convfctive  and  radiative  heating 

150  READ(5*3000)  (TITLE<L) *L=t ,12) 

IF(TITLE(1> .EO.PRVOUS)  GO  TO  POO 
READ ( 5*  3002 )  FCONV.FRAD 
WRITE <6* 155 )  FCONV,FRAD 

155  FORMAT ( 1H0  ,6HFC0NV=* 1PE12 .5  >4X5HFRADs» 1PF12.5/ ) 

PROPERTIES  OF  ABLATION  MATERIAL 
200  READ ( 5, 3000 )  ( TITLE (L ) *L=1 , 12 ) 

IF(TITLE(1). EO.PRVOUS)  GO  TO  300 
RFAO(5#3000)  (HEADNG(K) .Krl.lP) 

RFAD ( 5# 3002 )  TABL»TCHAR,TBEC.PHOV*RHOC»FRLOW,FMV,EMC,H30n.VL.HV, 

1 VPT , FV*TV ,CHARK , CHARC , ABLK , ABLC 
REA  (5*3003)  NP,NKC»NCPC.NKV.NCPV»NREC 
READ (5* 3002)  (TKC(K) ,XKC(K> ,Krl,NKC> 

RFAD ( 5  *  3002 )  (TCPC(M) , CPC (M) ,M=1 ,NCPC > 

RFAD ( 5* 3002)  (TKV(L) , XKV (I  ) *L=1 *NKV) 

RFAD(5*3002)  (TCPV(N) ,CPV(N) ,N=1.NCPV) 

RFAD(5*3002)  (TS(I) ,SR( I ) , 1=1 ,NRFC> 

WRITE (6*  3008 )  (HEADNg(K) ,K=1,12) 

WRITE (6* 210)  TABL»TChAR,TPEC,RHOV*RHOC,FRLOW,FMV,EMC,H30O,VL.HV* 

1 VPT  *  FV. TV  *CHARK *  CHARC  *  ABLK  *  ABt C 

210  FORMAT ( 1H0  *5HTABL=  »1PF12.5*3X»6HTCHAR=,1PF12,5*3X* 5HTREC=* 1 PE) 2.5* 
13X,5HRH0v=*1PF12.5*3x,5HRH0C=,1PF12.5,21X/1X,6HFbL0W=,1PE12.5,4X,4 
2HEMV=,lPei2.5»4X,4HEMC=.lPE12.5*3X*5HH300=*lPF12i5*5X,3HvL=*lPEl?. 
35/4X,3HHV=*1PF12.5»4x,4HVPT=.1PFi2.5,5X,3HFV=,1Pf12,5,5X,3HTV=,1PE 
ll2.5,PX»6HCHARK=,lPEl2.5/tX,6HCHARC=,lPEt2.5,3X,5HARLK=»lPE12.5>3X 
2 « 5HABLC=  * lPEt? . 5/ ) 

VLISVL 

VL=VL/12.0 

VLV=VL 

WRITE ( 6* 2P0 )  NP*NKC»NCPC,NKV,NCPV*NREC 

220  FORMAT ( 2X , 3HNP= . 1 14 . 4X »  4HNKC= , 1 1 4  *  4X , 5HN0PC= . 1 1 4 . 4X  *  4HNK V=  *  1 14. AX* 
15HNCPV= . 1 14  *4X . 5HNREC= . 114) 

WRITE (6.221 ) 

221  FORMAT(/32X,15HVIRGTn  MATFRI AL/20X.7HTHERMAL. 38X , 8HSPFCIFIC/3X . 1 1 H 
1 TEMPERATURE  *4X  * 12HC0N0UCT IVITY.19X. 11HTEMPERATURF.7X, AHHfAT) 

KLLL=MINO(NKv.NCPV) 

WRITE (6. 222)  ( TK V ( L ) , XKV (L ) .TCPV (L) *CPV(L) *L=1 .KLLL) 

222  FORMAT (2X »1PE12.5»4X,iPE12.5«18X, lPFlp .5  *3X,lPFl2.5) 

IF(NKV-NCPV)  223,227,225 

223  KLLLL=KLLL+1 

WRITE (6,224 )  (TCPV(L) , CPV (L ) .LsKLLLL, NCPv > 

224  FORMAT (48X»1PF12.5»3x» 1PF1 2 . 5 ) 

GO  TO  227 

225  KlLLLsKLLL+l 

WRITE (6, 226)  ( TKV ( L ) , XKV (L ) ,L=KLlLL*NK V ) 

226  FORMAT (2X»1PF12.5»4X,1PE12,5) 

227  WRITE(6.228) 

228  FORMAT ( //33X , 1 4HCHAR  MATFR I AL /2 OX , 7HTHERMAL , 38X , 8HSPFC IF IC/3X , 1 1 H 
1TFMPERATURE»UX» IPHCOnDUCTIVITY , lqX , 1 1HTEMPERATURf*7X , 4HHEAT) 

KLLL=MINO(NKC,NCPC) 

WRITE (6*222)  ( TKC (L ) , XKC ( L ) *TCPC ( L ) .CPC ( L ) *L=1*KlLL) 

IF(NKC-NCPC)  230,235,232 
230  KLLLL=KLLL+l 

WRITE (6,224 )  (TCPC(L) ,CPC (L ) ,L=KLLLL ,NCPC ) 

GO  TO  235 


r"  *  j.  au 

A1160 

A1170 

Alien 

A1190 

A1200 

A1210 

A1220 

A1230 

A1240 

A1250 

A1260 

A1270 

A1280 

A1290 

A1300 

A1310 

A1320 

A1330 

A1S40 

A1S50 

A1360 

A1370 

A1380 

A1390 

A1400 

A1410 

A1420 

A1430 

A1440 

A1450 

A1460 

A1470 

A148R 

A1490 

A1500 

A1510 

A1520 

A1530 

A1540 

A1550 

A1560 

A1570 

A 1580 

A1590 

A1600 

A1610 

A1620 

A1630 

A1640 

A1650 

A1660 

A1670 

A1680 

A1690 

A1700 
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232  Kt  LLU=KLLI  +1  A1710 

wPlT£<6,226)  <TKC<L) ,XKC(I  >  ,L=KL|  LL.NKC)  A17?0 

235  WRITE (6,240 )  A1730 

240  F  ORMAT  (  //?8X  •  P3HS1  IRF  ACE  RFCFSSION  T ABLE//R5V ,  1  IHTFMPERATiiRE ,  AY ,  1 1  H  A1740 
IBP  -  TN/SFC)  A1750 

WPITE ( 6  >245 )  (TS(T> ,SR<T) ,I=1,NRFC)  A1760 

245  F0RMAT(24X,1PF12.5»7x,1PF12.5)  A1770 

A17A0 

PROPERTIED  OF  TRAJECTORY  A17qn 

300  RFAO  ( 5*  3t)00 )  (  TIT LE  (L 1  * L  =  1  . 1?  )  AlftOO 

TF  (TITLE  ( 1  )  .FO.PRVO(IB)  Co  TO  400  Al«10 

RF AD ( 5 »  300 0  )  (HFAnNfi(L) ,t  =1*1?)  A1A?0 

RFAO ( 5  *  3004 )  NTRAPT  A1R30 

RF  AD ( 5 »  30 1 0 )  (TIMF(K) ,QCON(K ) ,QRaD(K) ,VFL (K) ,K=1 .NTRAPT)  Alfl40 

WRITE ( 6 »  3008 )  (HEAONr(L) ,1=1 ,12)  A1850 

wPlT£(6*310)  NTRAPT  A1P60 

310  FORMAT< 1 H0*27M  NO.  OF  TRAJECTORY  POINTS  =,114)  A1870 

wP ITE ( 6 ,  320  )  A1880 

320  F0RMAT(//PX,4WTIMF,«X, l?HO  CONVEcTI VE » 4X » 1 1HO  RAnTATl\/E,7X» BHVEL  OC  A1890 
1  IT  Y  )  AlO'OO 

WRITE <6, 330)  (TIMF(K) ,QCON(K) ,QRaO(K) ,VFl(K) ,k  =  1 .NTRAPT)  A 1910 

330  FORMAT  (1P4E1(S. 5)  A  19?0 

C  A1930 

C  PROPERTIES  OF  BACK-UP  STRUCTURE  A1940 

400  READ(5,3tlOO)  (TITLE(L)  ,L=1  ,121  A1950 

IF  (TITLF  (1  )  .FO.PKVOIIs)  Go  TO  500  A1960 

WRITE<6,410)  A1970 

410  format(//iox,31h  properties  of  backup  structure/)  ai9bo 

RFAD(5,3u07)  MMP,NPRs,8L  A1990 

RFAD(5,3002)  (XNPM(K) ,K=1 ,NMR)  A2000 

RFAD ( 5 , 415 )  (NKPR(I).MCPR<I),T=1,MMB)  A2010 

415  F ORMAT ( 1 U T5 )  A2020 

DO  420  K=1 ,NmO  A2030 

NPM(K)=YNPM(K)+0.000nO0n?  A2040 

420  CONTINUE  A2050 

wPITE ( 6 , 42b )  NMP,NPPS,BI  -  N  :  A20ft0 

425  F0RMaT(/4X,35HN0.  OF  MATERIALS  IN  BACK-UP  SHI FLD= , 1 I4/4X , 40HTOT A|  A207n 

INUMBEP  OF  NODES  in  Back-up  SHTEl  n=»lI4/4X*2BHTHICKNFSS  OF  BACK-UP  A2080 
2SHIFLn=,lPEl?.5//)  A2090 

RL=BL/12.0  A2100 

no  440  1=1 »nmr  A21 in 

LK=NKPB(1)  A2120 

ICP=NCPB(T)  A2130 

RF  AD ( 5 »  3u  00 )  ( (X1DNT(K,I) )*K=1  ,1?)  A2140 

READ  (  5, 3ii02 )  (  (TXK(J.I)  ,XK(J,T)  )  ,J=1,LK)  A2150 

READ ( 5, 3002 )  ( ( TCP ( J , I ) , CPX < J , I ) ) , J=1 ,LCP)  A2160 

WRITE ( 6 »  432 )  (XIDNT(K.I) ,K=1,12)  A2170 

432  fORMAT(//12Aa)  A21BO 

WRITE ( 6*  433 )  A2190 

433  FORMAT ( //?0X, 7HTHFRM aL » 38V » BHSPEc IFTC /3X, 11HTFMPfRATurE»4V,  1  ?HCOND  A220n 

l(iCTlVITY,19X,1  IHTFMPfrATi  IRE  »7X ,4hHE AT )  A2210 

KLLL=MINo(LK,l  CP)  A2220 

DO  434  N=) ,KLI  L  A2230 

WR 1 TE ( 6  *  222 )  (TXK (N, I ) ,XK (N, 1 1 .TcPIN, I)  ,cPX(N, I) )  A2240 

434  CONTINUE  A2250 

IF(LK-LCP)  435,440*4^7  A2260 

435  kLLLL=KLLI+1  A2270 


no  43ft  N=KLLU#ICP  A2280 

wP ITE ( 6 ,  224 )  <TCP(N,t) ,CPX(N,T) )  A2290 

436  CONTINUE  A2300 

CO  TO  440  A2310 

437  kLLLL=KLLL+l  A2320 

no  438  NzKLLLI  * LK  A2330 

WRI TE ( 6  *  226 )  ( TXK ( N , r ) , XK ( N *  I ) >  A2340 

430  CONTINUE  A2350 

440  CONTINUE  A236D 

RF  AD ( ft  *  3o02 )  (RHORX(L) *XRM(L> ,EMFP(L> ,EMrR(L) ,L=1 ,NMR)  A2370 

RF  AD ( ft  *  3002 )  (H( J) ,CAPX< J) ,FTFftT< J) ,BTFftT( J1 , Jrl.NMR)  A2380 

wP I TE  <  6  * 4ftO )  A2390 

4ftO  FORMAT (///55X* 10HFMI sft IVTTY/8X  » 8MMATEP I AL  *  5X «  7HDFNSI TY  *7X*9HTHICKN  A24n0 

lFftft*7X»  ftHFRONT  *  9x  * 4HRACK  > 7X  > 14HN0DES/MATFP I AL/ )  A24I0 

DO  460  LLJSl.NMP  A2420 

WRITE ( 6 »  486 )  LLJ,PHOnv(l L.J) .XPM(LLJ) > FMFp ( LL J ) *FmRB(LLJ> .XNPM(LLJ)  A2430 
4R5  FORMAT (11X.1I1 * 8X * lPF 1 0 . 4 , 4X . 1 PEl 0 . 4 , 4 X  *  1  PE 1 0 . 4 . 4X , 1PF1 0 . 4 , ftX , 1PF1  A244  0 

10.4/)  A2450 

4ftO  CONTINUF  A2460 

WRITE (6  * 4ft5 )  A2470 

4A5  FORMAT ( //4X  *  ftoHADDI T I ON A I  DATA  FOR  INDIVIDUAL  MATFRIAlS  IN  RACKtIP  A2480 
1 ftTRUCTURE //I l X  *  RhMATfR I Al  ,5X,16HpILM  COFFFIC IFNT , ftX , 13HCAP  THTCKmF  A2490 
2ftft  *  8X  * 5HE TEST.13X* 5HRTEftT )  A2500 

DO  400  Jr 1 , Nm°  .  A2ftl 0 

WRITE ( 6  * 470 )  J*  H(J) ,GAPX< J) ,FTFST( J) ,RTFST( J)  A2820 

470  FORMAT (13X*1tx*12X*1PF10.u»9X* 1PF 10.4*7X.t  PF 1 1 .4,7X, 1PE1 1 . 4/ )  A2B30 

480  CONTINUF  A2B40 

C  A2ftfi0 

C  PROPERTIEft  OF  ENVIRONMENT  A2B60 

600  RFAD (8*3000)  (TITLE(L) *L=1 *12)  A2570 

IF(TITLE(1 ) .FO.PRVOUR)  CO  TO  ftOO  A2580 

RF  AD  (  ft  *  3l(00  )  (HEADNGfL) *L=1*12>  A2ft90 

RFAD ( ft «  30  02  )  TENV.HFnV.FfNV.OI  Oftft  A2600 

WRITE (6  *  3008 )  (HFADNh(L) ,L  =  1.12)  A2610 

wP I  TE  ( 6  *  620 )  TENV  *  HFn V i FFNV  > QL  OSft  A2620 

520  FORMAT ( /4X  « 1 RHTFMPER ATUPF- *1PF12.8*4X, 17HFIL.M  COfFFIC TENT= . I PF 12. ft  A2630 
1.4X.12HVIFW  FACT0R=»1PE12.5*4X,7H0  L05T=. 1  PEI 2 . ft )  A2640 

C  A2650 

C  INITIAL  TFMPEPATuPE  niSTRTBUTTON  A2660 

600  RFAD ( ft »  300  0 )  (TITLE(L> *1=1 *12)  A2670 

IF(TITLE( 1 ) .FO.PRVOUft)  CO  TO  700  A2680 

RF  AD (ft.3000)  (HFAnNG(L) ,L=1.12)  A2690 

nPF— NP+NPPS  A2700 

TLzVL+BL  A2710 

XNPrNP  A2720 

OXrVL/ ( XnP—1 . 0 )  A2730 

nXX=DX  A2740 

RF  AD ( ft  *  30  02 )  TE8T?*TfmPi ,TX0  A27ft0 

IF ( TE8T2 )  61o, 620. 62o  A2760 

610  RFAD  (  ft  *  3(102  )  (TEMPI («) ,Kz1 *NPF)  A2770 

DO  61ft  K=1*NPF  A2780 

TX1(K)=TEMDI(K)  A2790 

TX2(K)=TX1 (K)  A280O 

TtJLl(K)  =  TXl(K)  A2810 

TIIL2(K)=TX1<K)  A2820 

bl 6  CONTINUE  A2830 

L=NP+1  A2840 
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no  6ip  1=1 *nmr 

I  N=NPM(I ) 

no  6i7  j=i , ln 

TX2T(J,I >=TFMnl(L) 

I  =L  +  1 

617  CONTINUF 

619  CONTINUE 
00  TO  62b 

620  CALL  TEMPO 

625  WRITE ( 6 » 30  08 )  < He ADNf,  <  L)  , I  =1 . 1 2  ) 

TF<TEST2)  630,636,6Un 
630  wP ITE ( 6 ,632 ) 

632  F0RMAT(4X,52HTEMPFRATURF  nlSTPIPllTlON  IN  HEAT  SHlFLD  TS  APBPITARV/ 
1  ) 

wPlTE (6,633)  (TFmDI <K1 *K=1 ,NPF) 

633  FORMAT11PAE12.5) 

GO  TO  64b 

635  wP ITE ( 6 ,637 )  TEMPI 

637  FORMAT ( //4X » 64  HT  EMPFR  A  TURF  DISTRrFUTION  IN  HEAT  PHIFLO  K  UNIFORM 
1  AND  EOUAL  TO  ,1PE10.U/) 

GO  TO  64b 

640  wP I TE ( 6 , 64 1 ) 

641  FORMAT <4X,54HI  INEAR  TFMPfRATUPE  OISTRIBUTTON  ASSUMED  TN  HFAT  SHlFL 
in/) 

WRITE ( 6 ,633 )  (TEMDIIL) ,L=t.NPF) 

645  IF(DMP)  700,700,646 

646  WRITE ( 6 , 647 ) 

647  FORMAT ( // ) 

640  wP ITE ( 6 , 649 )  (TXl (L) ,TX?(| ) ,L=1 ,NPF) 

649  FORMAT  (  2X  ,  1PF  1 2 « 5  »4X , 1 PFl 2,5 ) 

WRITE ( 6 ,650 ) 

650  FORMAT {//) 


enthalpy  as  a  function  of  tfmperatupe 

700  PFAD  (  6, 3(100  )  (TITLE(L)  ,L=1 ,12) 
TF(TITLE(1) .EO.PRVOUS)  GO  TO  726 
RFAD ( 5 »  3004 )  NHP 

RFAD l 6 , 3002 )  (HX(K) ,TW(K) ,K=l,NHp) 

725  DO  720  1=1 , Np 
TP(I)=0 
tpi ( I i =o 

TR2 ( I ) =0 
TFM ( I ) =0 
XMDG ( T ) =0 , 0 
728  CONTINUF 

WRITE ( 6 , 730 ) 

730  FORMAT ( 1H1 ,12H0UTPUT  nATA.//) 

XC( 1)=0.0 
no  740  1=2, NP 
XC(I)=XC(T-1)+DX 
740  CONTINUF 

760  IFIT-TIME(NK) )  766,770,760 
760  NK=NK+1 

TF(NK-NTRAPT)  750,750,762 

762  WRITE ( 6 , 763 )  MK 

763  FORMAT ( 1HO , 33H  THE  VALUF  OF  NK 
GO  TO  905 


A2R50 
A206O 
A2870 
A2A80 
A209O 
A2900 
A2910 
A2920 
A2930 
A2940 
A2950 
A2960 
A2970 
A2980 
A2990 
A3000 
A3010 
A3020 
A3030 
A3040 
A3050 
A3060 
A3070 
A3O0O 
A3090 
A3100 
A31  10 
A3120 
A3130 
A3140 
A3150 
A3160 
A3170 
A3180 
A3190 
A3200 
A3210 
A3220 
A3230 
A3240 
A3P50 
A3260 
A3270 
A3280 
A3290 
A3300 
A  331 0 
A3320 
A3330 
A3340 
A3350 
A3360 
A3370 
A  3300 
A3390 

IS  IN  FRROP,  NK= , 1 14 )  A3400 

A3410 
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765  IF ( NK-2 )  762,766.766  A3420 

766  0C0NX=QC0N(Nk-1>+(  (OcON(NK)-OOON(NK-l)  > / < TIME ( NK ) -TIMF < Nk-1 )  )  )  A3430 

1*(T-TTMF(NK-1 ) )  A3440 

oconx=fconv*qconx  A3450 

OPADXrQRAn(NK-l)+( (ORAD(NK)-OPAn(NK-l ) > / ( TIME ( NK ) -TIMF ( Nk-1 ) ) )  A3460 

1* ( T-TTMF ( NK-1 ) )  A 3470 

(5RADX=FRA0*0RA0X  A34A0 

VFLX=VEL(NK-1 )+( (VEL(NK)-VEL(NK-l 1 )/(TIMf(NK)-TImF(NK-1> ) )  A3490 

1* ( T-TTMF (NK-1 ) )  A3500 

GO  TO  775  A3510 

770  OC0NX=FC0NV*0C0N(NK)  A3R20 

ORA0X=FRAn*QRAD(NK )  A3530 

VFLX=VEL(NK)  A3540 

C  A3550 

C  COMPUTE  HFAT  PLOCKAGf  AT  FRONT  SllPFACF  A3560 

775  TF ( I 17-1 )  77A.77fl.776  A3570 

776  TF ( I 17-NHP )  777,777,778  A35A0 

777  TF(TX?(INT)-TW(I17) )  782.78fl.7flO  A3590 

778  WRITE ( 6 »  779 )  TX2(INT)  A3600 

779  FORMATdHn.floM  THF  RAMGF  OF  THE  FNTHAL PY-TEMPFR AtURE  CURVF  FIT  WAR  A3610 

IfXCEEDEO  AT  A  TEMPERaTURF  0F,lEln.4)  A3620 

GO  TO  906  A3630 

7A0  tl 7=117+1  A3640 

GO  TO  77h  A 3650 

7fl2  TF(TXP(I NT )-TW (117-1) 1  7fla.7flfl.7fl6  A3660 

7P4  117=117-1  A3670 

GO  TO  776  A3680 

7P6  hW=HX< 117-1 )  +  ( (HX(I17)-HX(I17-1) > / ( TW ( 1 1 7 ) -TW ( 1 17-1 ) ) )  A3690 

1*<TX2(INT )— T  w ( 1 17-1 ) )  A3700 

GO  TO  7A9  A3710 

7fl8  HWzHX (117)  A3720 

7A9  mTX=H300+( (VF!  X**2 ) /SO056 .5 )  A3730 

OPLOCK=(FRLOw*XMOG(lNT)*(HTX-HWn/3600.0  A3740 

C  A3750 

C  COMPUTE  HFAT  IN  DUE  TO  SllPFACF  COMBUSTION  A3760 

XMOO=XMDC  A3770 

CALL  OXIUAT(XMDO.OOXrn)  A37AO 

C  A3790 

C  COMPUTE  Q-HOT  WALL  A3AOO 

XFUOMP.EO.O.  )  GO  TO  U001  A3B1 0 

TF(T.GE.TDMP)  DMprl.n  A3A20 

4001  7- ( HTV-Hw ) / ( HTX-H300 )  A3A30 

TF(Z-l.O)  79n,792, 793  A3A40 

790  I F  (  Z )  791,791,793  A3ARO 

791  OHW=0.0  A3A60 

GO  TO  1790  A3A70 

792  f)HW=QCONX  A3P80 

GO  TO  1790  A3A90 

793  OHW=Z*QCONX  A3900 

1790  7ZZ=(0HW-0BL0CK)/0HW  A3910 

TF(ZZ7-0.2)  1798, 179fl, 179U  A39?0 

1798  oPLOCK=0 . a*Ohw  A3930 

C  A3940 

C  NFT  HPAT  INTO  FRONT  SURFACE  A3950 

1794  IF ( IEM { InT ) )  795,795,797  A3960 

795  TF ( TX? ( INT )— TCHAR )  706,796,797  A3970 

796  FMX=EMV  A390O 
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GO  TO  798  #3990 

797  IFM(INT)=1  A4000 

FMX=EMC  A4010 

798  GlN=0RA0x+QHw+Q0XI0-ORL0CK-<4.8333E-13)*FMX*FV*( (TX2 { TNT ) **4 )-  A4020 

1 ( TV**4 ) )  A4030 

IF(OMP)  804#8o4»800  A40t*O 

800  wPITE (6»80l )  A4050 

801  FORMAT ( /// )  A4060 

WRITE (6»B02)  OCONX , OR ADX , VELX  »HTX  »HW»  ? ,GRLOCK  , QHw *QOXTD,GTN  A4070 

802  F0RMAT<1X,6HQ00NX=,1PF17.5.2X,6H0PADX=,1PF12.S»2X.5HVFLX=,1PF12.B,  A4080 

l?X,4HHTX=,lPEl2.5,2X,3HHw=.lPF12.5/lX,2H7=*tPF12.5.2X,7HQRL0CK=»lP  A4nq0 
2f12.5,2X,4HQHW=,1PE12.5,7X,6H00Xi0=»1PE17.5,2X,4h0IN=,1PF12.5/)  A4100 

804  OTN=QTN*3600.0  A4110 

A4120 

CHECK  FOR  PROMT  SURFACE  RFCESSION  (CHAR  LAYFR  RFMOVAL)  A4130 

CALL  RECESS ( XMDC  , XLOsT  »TRFC»DT»PH0C,TS»SR»TX2(1)»  NRFC  , NRS, ERRS  »SvO  A4140 
1 • SCOT , OMP )  A4150 

IF ( ERR5 )  8050, 8060. 9nS  A4160 

8060  Vl  V=VLV-XI  OST  A4170 

XI  STV=XLSTV  +  xt  OST  A4180 

XI  STI=XLSTV*12.0  A4190 

nXV=VLV/(XNP-1 .0)  A4200 

XV< 11=0.0  A421 0 

no  1780  1=2, NP  A4220 

XV(I)=XV(T-l)+DXv  A4230 

1780  CONTINUE  A4240 

nX=DXV  A4250 

IFIFRR4)  806, 806. 80S  A4260 

80b  00  TO  90b  A4270 

806  CALL  COFFF (NpFT ,SOOT )  A4280 

IF(OMP)  8069,8069,8001  A4290 

8061  WRITE ( 6  #  8062 )  A4300 

8062  FORMAT ( /IX , 23H  COEFFICIENTS  FOR  SWUFT/)  A4310 

nO  8066  I=1*NPFT  A4320 

WRITE ( 6 , 8064 )  A ( I ) , 8 ( I ) , C ( 1 1 , D ( I) , I  A4330 

8064  FORMAT (1H0,5HA(I)=»1 PF12 . S , 2X , 5Hp ( 1 1  = , 1PF1 2  .  S , 2X , SHC  C  T )  —  . 1 PF 12,8,2  A43U0 

1X,5HD(I)=,1PF12.5,2X,2HI=,I3>  A43B0 

8066  CONTINUE  A4360 

8069  IF ( FRR2  )  807,807,808  A4370 

807  TF(ERR3)  810,810.808  A4380 

808  wPITE ( 6 » 809 )  TKK  A4390 

809  FORMAT (1H0.18H  THF  VaLUF  OF  IKK=,1I4)  A4400 

60  TO  90b  A4410 

810  CALL  SWUFT (A,R,C,D»Ty, NPFT , DMP )  A4420 

827  no  828  1=1, NP  A4430 

TX1 ( I ) =TX2 ( I )  A4440 

TX2(I)=TY(I)  A4450 

828  CONTINUE  A4460 

CALL  nON2 (XLOST *XV»TX2*  Np ,XC ,TX2C  *  XDV ,KK\/,XLSTV»nXX)  A4470 

830  CALL  ABLATE  A4480 

XMDT=XMnG ( INTI +XMOC  A4490 

LT=NP+1  A4S0O 

no  1815  1=1, NMB  A4S10 

LLT=NPM(1>  A4520 

TF(I.FQ.l)  GO  TO  1812  A4530 

TF(GAPX( 1-1) .FQ.O. )  GO  TO  1812  A4540 

KKT=1  A4S50 
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GO  TO  1*13  A4560 

1812  KKT=2  A4570 

1813  00  1815  JrKKT.LLT  A4580 

TX2T  ( J,  I )  =TY  ( LT  )  A4590 

I T=LT+1  A4600 

1815  CONTINUF  A4610 

no  1819  I=1»NMB  A46P0 

TF(I.FO.l)  60  TO  1  SI  ft  A4630 

TF<GAPX(I-l).FQ.O.  >  GO  TO  1817  A4640 

GO  TO  1819  A4650 

1816  TX2T  { 1  #  1  )  =TY  ( KIP)  A4660 

GO  TO  1819  A4670 

1817  I  X=NPM(I-1 )  A4680 

TX2T( 1 , T )=TX2T(lX, 1-1 )  A4690 

1819  CONTINUF  A4700 

I M=NP+1  A4710 

no  833  t=1»NMP  A4720 

|.7=NPM(I)  A473D 

no  833  J=1 » L?  A4740 

TV2(LM)=TX2T< J.I)  A4750 

I  M=LM  +  1  A4760 

833  CONTINUE  A4770 

00  5834  1=2 . NPTT  A4780 

TF ( T-TT ABLE (II)  5835,5835,5834  A4790 

5835  nTS=DFLTT  ( 1-1  )  A48nfl 

TPRCT=IPHC (1-1 )  A4810 

nT=OELTT(T-l 1/3600.0  A4820 

GO  TO  5836  A4830 

5834  CONTINUE  A4840 

nT5=UPLT1 (NPTT)  A4850 

t  PRCT=IPKC ( NPTT )  A4860 

f1T=0El  TT ( NPTT 1/3600,0  A4870 

5836  tCT=ICT+l  A4880 

5838  vl  TEM=SAVY3  A4890 

CALL  TSOTHM<XV.TX2.1060,  ,NP ,5 AVEIT)  A490O 

9AVFIT=5AVEIT+XLSTV  A4910 

TF( SAVY3.LT. SAVF IT) 5a VY3=5AVETT  A4920 

TF(VLTEM.EQ.5AVY3)G0  TO  839  A49J0 

5AVX=T  A4940 

5A V YlrXLST I  A4950 

5AVY2=TX2 (NP)  A4960 

CALL  TSOTHMIXV.TX?. 1460. .NP.SAVYu)  A4970 

839  Rl  TFM=SAVY4X  A49B0 

CALL  TSOTHM ( X V , TX2  f 1 460 • » Np .WFKEpP)  A4990 

wEkFEP=WEKEEP+XLSTV  A5000 

TF(  5AVY4X.lt.  wEKEEPKaVY4X=WEKEEP  A501  0 

TF(BLTEM.EQ.5AVY4X)G0  TO  838  A5020 

5AVEXX=T  A5030 

6AVY1X=XL5TI  A5040 

9AVY2X=TX?(NP1  A5050 

CALL  TSOTHM (XV»TX2»l060.« NP , 5 A VY3X )  A5060 

838  CONTINUE  A5070 

IF(IPRCT-TCT)  835,835,840  A5080 

835  WPITE (6 , 837  >  T , GCONX , oRAnx , VEI  X , XMDG( TNT) , XMDC , XmOT , XI  ST I , QHW  A5090 

837  EORMAT(1HO,5HTIME=,  A5100 

1  1PE1 2.5.2X , 12HOCONVECTT VE= , 1  PE 1?.5,?X,11 HOR ADT  AT  A5110 

1 T VE- , 1  PEI? . 5, ?X,9MVFlOCITY=,1PE12.5/1X,1 flHGAS  APlATIOM  RATEr.lPEl?  A5120 


2 . 5 » 2X  » 19HCHAR  ABLATION  P ATE=, 1 PFl 2 . 5 . 2X , 20HTOTAL  ABLATION  RATF=»1P  A5130 
3F12.5/1X. 16HrfCFSSI0n  OFPTH= , 1 PFl 2 . 5 t 2X » 1 OHOHOT  wALL= , lPFl 2. 5)  A5140 

840  T=T+OTS  A5150 

841  TF(NPLOT.NE.l )  GO  TO  842  A5160 

CALL  SAVE ( ASAvEl , ASAvF2» ASAVF3.UsFAf NX A , XL  ST I , DTs t TLIM.T , VALIIFA )  A5170 

CALL  SAVE(BSAVEl,RSAvF2,BSAVF3.UsFB.NXB,TX2(NP) ,DTS*Tl IM.TrVALUFR)  A51B0 

CALL  IS0THM(XV*TX2f  1(160.  ,NP,Y3)  A5190 

CALL  SAVE ( CSAVE1 »CSAvF2»CSAVF3*UsFC»NXC»Y3»DTS.T|_IM»T» VALUEC )  A5200 

CALL  ISOTHM(XV»TX2*l460. .NP,Y4)  A5210 

CALL  SAVE ( DSAVE1 » OSA VF2. HS AVF3.USFD»NXD»Y4, DTS .TLTMiT, VALUED)  A5220 

T F  ( USF A « NF • 0 . 0 ) GO  TO  9842  A5230 

I F ( USFB , NF . 0 . 0 ) GO  TO  9842  A5240 

TF(USFC,NF.0.n)GO  TO  9842  A5250 

T F ( USFD . NF • 0 , 0 ) GO  TO  9842  A5260 

GO  TO  9843  A5270 

9842  xPLOT=T-DTS  A5280 

YPLOTt=VALUEA  A5290 

IF  ( USFA.  NF  •  0  .  (1 )  YPLOT 1  zU^F A  A5300 

YPLOT2=VAl UER  A5310 

IF (USFB .NF , 0 . n ) YPLOT  2FUSFP  A5320 

YF’L0T3=VALUEC  A5330 

T  F ( USFC . NF . 0 . 0 ) YPLOT^rUSEC  A5340 

YPLOT  4=V AlUFn  A5350 

IF(USFD.NF.O.O) YPLOTarUSFn  A5360 

WRITE  (ll)XPLOT»YPLOT1 .YPI OT2.YPI 0T3.YPL0T4  A5370 

9843  IF(ICTP.NF.O)  GO  TO  842  A5380 

ICTP=1  A5390 

XPLOTrT  A5400 

YPL0T1=XLSTI  A5410 

YPL.0T2=TX2(NP1  A5420 

CALL  TSOTHM(XV.TX2. 1060. ,NP»YPL0t3>  A5430 

CALL  TSOThM(xV,Tx2. 1460. .MP.YPL0T4)  A5440 

WRITE  (11) XPlOT  t YPLOT 1 • YPLOT2  * YPLOT 3 »  YPLOT4  A5450 

842  IF ( IPPCT-ICT  )  845,845,900  A5460 

845  wPITE ( 6  , 850  )  T  A5470 

TPCT=TPCT+1  A5480 

lF(IPCT.EO.2)TPCT=0  A5490 

lF(IPCT.EO.0)TCTP=0  A5500 

8so  format<iho,7phtfmperaturf  distribution  in  hfat  shield  at  thf  fnp  o  assio 
IF  THE  TIMF  STFP,  T=  , 1  PEI  2 .5 , 1 X , 7HSFCONDS// >  A5520 

WRITE  <  6 « 860 )  A5530 

880  FORMAT ( 4X  » 49HTEMPFR A  TURF  DISTRIBUTION  IN  THF  ABLATING  MATFRI Al //  I  A5540 
kKV=KKV+1  A5550 

WRITE ( 6»  882 )  (TX2CII) .1=1 *KKV)  A5560 

882  FORMAT ( 6X  * 1PF 12.5, IPSF 18 . S )  A5570 

IJ=NP+1  A5580 

WRITE (6  *  884 )  A5590 

884  F0RMAT(//4X.aoHTEMPFRATURF  DISTRIBUTION  TN  THP  BACK-UP  STRUCTURF//  A5600 
1)  A5610 

WRITE  ( 6  •  882 )  (TX2(I)  ,I  =  IJ,NPF)  A5620 

WRITE ( 6 » 885 )  A 5630 

865  FORMAT!//)  A5640 

ICT=0  A5650 

900  CONTINUE  A5660 

IF (T-TLIM)  760  ,750  » 90S  A5670 

905  TF(NPLOT.NE.l)  GO  TO  909  A5680 

XAVY3=SAVY3-SAVYl/12.  A5690 
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XAVY4X=SAVY4X-SAVY1X/12.  A5700 

TF(SAVX.EQ.XPL0T)G0  TO  9005  A5710 

wPITE  ( 11 ) SAVX . SAVY1 *  SAVY2 ,XAVY3 *  9AVY4  A5720 

9005  TF(SAVEXX.E0.XPL0T)G0  TO  9006  A5730 

9AV4I=SAVY4X*12.  A5760 

9006  9AV3I=SAVY3*12.  A5750 

WRITE ( 1 1 ) SAVEXX*SAVY IX*  SAVY2X *  SAVY3X *  XAVY4X  A5740 

WRITE (6*929)SAV3I *SA V4I  A5770 

929  FORMAT ( 1H0  *23HMAXIMUM  1060  ISOTHfPM  =F16.8»?X23HmAXIMUM  1460  TSOTH  A5780 

IfRM  =F16.8)  A5790 

WRITE  (11) STOP*  STOP  * 9T0P*ST0P  r STOP  A5800 

909  IF (LPLOT .NE.NCASEIGO  TO  911  A5810 

DATA  FND/6H  FND  /  A5820 

wP I TE ( 1 1 ) FND , FND , END , FNn . FND , FND , FND*  FND . FND , FND . FND , FND  A5830 

OIIIT=8888.  A58U0 

WRITE ( 11 ) QUIT, QUIT *Ql)IT*OUIT,OUlT  A5850 

FND  FILE  11  A5860 

RFWIND  11  A5870 

911  TF(TE<;T2)  910*930*930  A5880 

910  DO  920  JJK=1,NPF  A5890 

TX1 (JJK)=TEMDT < JJK)  A5900 

TX2(JJK)=TX1(  JJK)  A  59.10 

TUL1(K)=TX1(K)  A59P0 

TUL2(K)=TX1(K)  A5930 

920  CONTINUE  A5940 

TL=NP+1  A5950 

DO  926  I=1*NMR  A5960 

TIN=NPM(I)  A5970 

DO  924  J=1*ILN  A5980 

TX2T(J.I)=TFMDI(IL)  A5990 

TL=1L+1  A6000 

924  CONTINUE  A6010 

926  CONTINUE  A6020 

GO  TO  940  A6030 

930  CALL  TEMPO  A6040 

940  t=TINT  A6050 

nX=OXX  A6060 

DTS=DFLTT(1)  A6070 

nTsDELTT  ( 1 )  /.3600. 0  A60B0 

VI  V=VL  A6090 

GO  TO  50  A6100 

END  A6110 
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SI8FTC  CO£F  B0000 

C  THIS  SUBROUTINE  DETERMINES  THE  COEFFICIENTS  OF  THE  MATRIX  B0010 

SUBROUTINE  COEFF ( NPFT > SDOT )  B0020 

C  B0030 

DIMENSION  TITLE  ( 12)  »HEADNG(12)  , XIDNT  ( 12 1 12)  rTKC  (20 )  t  XKC  <  20 )  #  B0040 

ICPC (20 ) *TKV (20 ) *  XKV (20) »TCPV(20) #CPV(20) # TIME (300) *QCON(300) »  B0050 

2QRAD ( 300 ) » VEL ( 300 ) » XNPM ( 12) »NKPB ( 12 ) # NCPB ( 12) >TXK(20,12) ,XK(20,12)  B0060 

3»TCP(20f 12) fCPX(20f 12) ,RHOBX ( 12) * XBM ( 12 ) » EMFB ( 12 ) »EMBB(12) *HXX(12)  B0070 

4  »GAPX (12) #FTEST(12J rBTEST ( 12) > TEMDI (200 ) * TX1 ( 200 ) # TX2 (200 ) »  B0080 

5TX2T ( 10 » 12) »TUL1 (200 ) »TUL2 (200 ) »HX (50 ) »Tw (50) • IR (50 ) ►  IR1 (50) •  B0090 

6IR2 ( 50 ) *  TUL ( 50 ) »IEM(50) > TY (200 ) » A ( 200 ) »B ( 200 ) *C (200 ) ,D ( 200 ) »  B0100 

7R ( 50 ) t RHO ( 50 ) »CP(50) »DXB(12) »XKB(10»12) >CPB ( 10 » 12) » XMDG ( 50 ) »  B0110 

8YK(50) » AB(10»12) »BB(10#12) »CB(10*12) *DB(10»12) » SB (10, 12) ,  B0120 

9RB1(10,12) *  RB2 (10,12) ,H(12) » S ( 50 ) »NPM ( 12 )  B0130 

DIMENSION  TTUL (50 ) ,RH0Y1 (50 ) »RHOY2 ( 50 )  » DRHO (50 ) ,  TCPC (20 )  B0140 

C  B0150 

COMMON  TKC»XKC> TCPC, CPC, TKV»XKV,TCPV,CPV,XNPM,RHOBX, XBM, EMBB,  B0160 

1EMFB,NKPB,NCPB,TXK,XK,TCP»CPX,NPM,GAPX,FTEST»BTEST,TEMDI,TX1,  B0170 

2TX2,TX2T  , TUL, TUL1 *  TUL2, IR » IR1 , IR2 » A»B»C  »D»S,R, AB,BB,CB,DB,SB,  B0180 

3RB1 , RB2 , T Y , RHOY1 , RH0Y2 ,  XMDS , RHO » CP » YK , X«B , CPB * DXB , DT , XLOST ,  B0190 

4TABL,TCHAR,TREC»RHOV,RHOC*FBLOW»EMV,EMC>H300»NKC,NCPC,NKV,NCPV,  B0200 

5NP » NMB » NPBS » NPF » TEST2 ,  TEMP I » TXO » TENV » HENV » FENV » QLOSS » TL I M  r  T I NT  B0210 

COMMON  I1,I2,I3,I4»I5»I6»QIN»INT »DX, XMT , TL» VL,BL ,DMP,ERR1 ,ERR2,  B0220 

1ERR3,ERR4,HV,VPT»CHARK,CHARC»ABLK»ABLC»XMDC,H  B0-230 

C  B0240 

CALL  PROP  B0250 

YNP=NP  B0260 

S(INT)=(RH0(INT)*DX*CP(INT> )/(2.0*DT)  B0270 

R(lNT)=(1.0)/( (DX/2.0) *( (1.0/YK(INT) ) + ( 1 . O/YK { INT+1 ) ) ) )  B0280 

A(INT)=0.0  B0290 

8 ( INT ) = ( - ( ( XMDG ( INT ) +XMDC ) *CP ( INT ) +S ( I NT ) +R ( INT ) -RHO ( INT ) *CP ( INT )  B0300 

1* (SDOT / (2 • 0* ( YNP-l . 0 ) ) ) ) )  B0310 

C ( INT ) =XMDG ( INT+1 ) *CP ( INT+1 ) +R ( INT) +RHO ( INT+1 ) *CP ( INT+1 ) *SDOT  B0320 

1*( (YNP-l. 5)/(YNP-1.0) )  B0330 

D ( INT ) = ( - ( QIN+S ( INT ) *TX2 ( INT ) ) )  + (XMDG ( INT) -XMDG ( INT+1 ) ) *HV  B0340 

NPP=NP-1  B0350 

JNT=INT+1  B0360 

DO  10  I=JNT,NPP  B0370 

XI=I  B0380 

S(I)=(RHO(I)*DX*CP(I))/DT  B0390 

R(I)=(1.0)/( (DX/(2.0*YK(I) > )+(DX/(2.0*YK(I+l) ) ) )  B0400 

A ( I ) =R ( I— 1 )  B0410 

B( !)=(-( XMDG ( I ) *CP ( I ) +R ( 1-1 ) +R ( I ) +S( I ) +  RHO ( I ) *CP ( I ) +SDOT* ( ( YNP-X I  B0420 

l-0.5)/(YNP-1.0)>)>  B0430 

C ( I ) =XMDG ( 1+1 ) *CP ( 1+1) +R ( I ) +RHO ( 1+1 ) *CP ( 1+1 ) +SDOT* ( ( YNP-XI-0 .5)  B0440 

1/ ( YNP-l • 0 ) )  B0450 

D(I)=(-(S(I)*TX2(I) ) )  + ( XMDG ( I ) -XMDG ( 1+1 ) ) *HV  B0460 

10  CONTINUE  B0470 

R(NP)=(1.0)/( (DXB(l)  /(2.0*XKB(1*1) ) )+(DXB<l)  / (2 . 0*XKB (2 , 1 ) ) ) )  B0480 

S ( NP )  =  < RHO ( NP ) *DX*CP ( NP ) +RHOBX ( 1 ) *  CPB (1,1) *DXB ( 1) ) / ( 2. 0*DT)  B0490 

A(NP)=R(NP-1)  B0500 

B(NP)=(-(XMDG(NP)*CP(NP)+R(NP-1)+R(NP)+S(NP) ) )  B0510 

C(NP)=R(NP)  B0520 

D (NP) = ( -S (NP ) *TX2 (NP) )  +XMDG(NP)*HV  B0530 

DO  200  1=1, NMB  B0540 

IF(I-l)  20,20, 30  B0550 

20  AB(1,I)=A(NP)  B0560 

BB(1»I)=B(NP)  B0570 

CB(1,I)=C(NP)  B0580 

DB(1»I)=D(NP)  B0590 
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GO  TO  65  B0600 

30  L=NPM(I“1>  B0610 

IF (FTEST ( I ) )  45*40*45  B0620 

40  SB(1*I)=(RH0BX(I>*CPB(1*I)*DX8(I)+RH0BX(I-1)*CPB(L*I-1>*DXB(I-1) )/  B0630 

1(2.0*DT)  B0640 

RB1(1*I)=(1.0)/((DXB(I“1)/(2.0*XKB(L*I“1)) >+(DXB(I-l)/(2.0*XKB(L-l  B0650 

1 r 1-1 ) ) ) >  B0660 

RB2(1*I>=(1.0)/((DXB(I)/(2.0*XKB(1*I) ) ) + (DXB ( I > / (2, 0*XKB (2* I ) )) )  B0670 

AB(1*I)=RB1(1*I)  B0680 

BB(1*I)=(“(RB 1(1*1) +RB2 (1*1) +SB (1*1)))  B0690 

CB ( 1  *  I ) =RB2 (1*1)  B0700 

DB(1*I)=(-(SB(1*I)*TX2T(1*I) ) )  B0710 

GO  TO  65  B0720 

45  IF (FTEST ( I )  )  50*40*55  B0730 

50  G=(1.73E-09)/(1.0/EMBB(I-l>+1.0/EMFB<I)“1.0)  B0740 

GO  TO  60  B0750 

55  G=0.0  B0760 

60  SB ( 1  *  I ) = (RHOBX (I)*CPB(1*I)*DXB(I) )/ (2 • 0*DT )  B0770 

RB2 ( 1  *  I )  =  ( 1 . 0 ) / ( (DXB(I)/(2.0*XKB(1*I)) >+<DXB(I)/(2.0*XKB<2*I> ) ) )  B0780 

AS ( 1  *  I ) =H ( 1-1 ) +4 . 0*G* ( TX2T <  L  *  1-1 ) **3 )  B0790 

BB(1*I)=(-(H(I-1)+4.0*G*(TX2T<1* I  )**3)+RB2(l»I)+SB(l*I> ) )  B0800 

CB ( 1  *  I )=RB2( 1»I)  B0810 

DB(1* I)=3.0*G*( (TX2T(L*I“1)**4)“(TX2T(1*I)**4) )-SB(l*I)*TX2T(l,I)  B0820 

65  LF=NPM(I)-1  B0830 

DO  100  J=2 * LF  B0840 

SB ( J* I ) =  (RHOBX (I)*CPB(JfI)*DXB(I)  )/DT  B0850 

RBl(J#I)=(1.0)/( (DXB(I)/(2.0*XKB(J-1»I)) ) + (DXB ( I ) / (2. 0*XKB ( J» I ) ) ) )  B0860 

RB2  («J»  I )  =  ( 1 . 0 )  /  (  (DXB ( I )  /  ( 2 •  0*XKB  ( J+l  r  I )  )  )  +  < DXB  ( I ) / (2. 0*XKB  ( J»  I )  )  )  )  B0870 

AB(JrI)=RBl(J'I)  B0880 

BB ( J* I )= (“(RBI ( I ) +RB2 ( J* I ) +SB( I ) ) )  B0890 

C8(J*I)=RB2(J,I)  B0900 

DB(J»I)  =  (-(SB( J»I)*TX2T(J»I) )  )  B0910 

100  CONTINUE  B0920 

IF(I-NMB)  110*250  *250  60930 

110  LNF=NPM(I)  B0940 

IF (BTEST ( I ) )  120.115*120  B0950 

115  SB ( LNF  *  I )  =  ( RHOBX ( I ) *CPB (LNF  *  I ) *DXB ( I ) +RHOBX ( 1+1 ) *CPB ( 1 . 1  +  1 ) *DXB ( 1+  B0960 

11>)/(2.0*DT)  B0970 

RBI ( LNF  *I)  =  (1.0)/( (DXB(I)/(2 . 0*XKB (LNF-1 *!)))  +  ( DXB(I) / (2. 0*XKB (LNF  B0980 

1*1))))  B0990 

RB2(LNF*I)=(1.0)/( (DXB (I+1)/(2.0*XKB( 1*1+1) ) ) + ( DXB ( 1+1 ) /  B1000 

1(2.0*XKB(2*I+1) ) ) )  B1010 

AB ( LNF  *  I ) =RB1 ( LNF  *  I )  B1020 

BB(LNF.I)=(“(RB1(LNF.I)+RB2(LNF.I)+SB(LNF*I) ))  B1030 

CB(LNF.I)=RB2(LNF*I)  B1040 

DB ( LNF • I )  =  ( - ( SB  <  LNF  *  I ) *TX2T ( LNF  *  I ) ) )  B1050 

GO  TO  200  B1060 

120  IF (BTEST ( I ) )  125.115*127  B1070 

125  G=  ( 1  »73E“09)  /  ( 1 . 0/EMBB  ( I )  +1 .0/EMFB  ( 1+1 )  “1  *0 )  B1080 

GO  TO  130  B1090 

127  G=0.0  B1100 

130  SB(LNF*I)=(RHOBX(I)*CPB(LNF*I)*DXB(I) >/(2.0*DT>  B1110 

RBI (LNF  *I)  =  (l*0)/( (DXB(I)/(2.0*XKB(t NF-l.I)) )  +  (DXB ( I ) / (2. 0*XKB ( LNF  B1120 

1*1))))  B1130 

AB(UNF.I)=RB1(LNF*I)  B1140 

BB (LNF .!)  =  (-( RBI (LNF* I ) +H ( I ) +SB (LNF *I)+4,0*G*(TX2T (LNF *  I ) **3) ) )  B1150 

CB ( LNF . I ) =H ( I ) +4 . 0*G* ( TX2T (1*I+1)**3>  B1160 

DB ( LNF  »I)=3.0*G*( ( TX2T  <1.I  +  1)**4)“(TX2T( LNF  *  I ) **4) )-SB ( LNF* I ) *TX2T  B1170 
l(LNF.I)  B1180 

200  CONTINUE  B1190 


250  MN=NPM(NMB)  B1200 

IF ( QLOSS )  270  r 260  *  270  B1210 

260  SB(MNrNMB>=(RHOBX(NMB)*CPB(MNrNMB)*DXE(NMB) )/(2.0*DT)  B1220 

RBl(MNrNMB)=(1.0)/( (DXB(NMB)/(2.*XKB(MN-1»NMB) ) )+(DXB(NMB)/(2.0*XK  B1230 

lB(MNrNMB)  )  ) )  B1240 

AB(MNrNMB)=RBl(MNrNMB)  B1250 

BB(MNrNMB)=(-<RBl(MNrNM8)+SB(MNrNMB> ) )  B1260 

CB(MNrNM6)=0.0  BI270 

DB (MN  r NMB ) = ( - ( SB ( MN  r NMB ) *TX2T ( MN  r NMB ) ) )  B1280 

GO  TO  280  B1290 

270  SB ( MN r NMb ) - ( RHOBX ( NMB ) *CPB ( MN  r NMB ) *DXB { NMB ) >/(2.0*DT)  B1300 

RBI (MNr  NMB)  =  ( 1 . 0 ) / ( (DXB (NMB) / (2. 0*XKB (MN-1 r NMB) ) )  +  (DXB (NMB) / (2 . 0*X  B1310 

1KB(MN»NMB) )))  B1320 

AB (MN»NMB)=RB1 (MNr NMB)  B1330 

BB ( MIM  r NMB )  =  (-( RBI (MNr NMB > +HENV+ < 1 . 73E-09) *FENV*4 . 0* < TX2T ( MN r NMB)**  B1340 

13)+Sb(MN»NMB) ) )  B1350 

CB  ( MN r NMB )  =0  •  0  B1360 

DB(MNrNMb)=(-(HENV*TENV+FENV*(1.73E-09)*( ( TENV**4 ) +3 . 0* < TX2T (MN r NM  B1370 

1B)**4) ) +SB ( MN  r NMB  >  *TX2T ( MN » NMB ) ) )  B1380 

260  L=NP+1  B1390 

DO  300  1  =  1  r NMB  B1400 

K=NPM(I)  B1410 

IF(I.EQ.l)  GO  TO  282  B14'20 

IF(GAPX(I-1) .EQ.O.)  GO  TO  282  '  B1430 

KT=1  B1440 

GO  TO  285  B1450 

262  KT=2  B1460 

265  DO  290  J=KTrK  B1470 

A (L ) =  AB ( J  r  I )  B1480 

B(L)=BB(JrI>  B1490 

C(L)=CB(J»I)  B1500 

D(L)=DB(JrI)  B1510 

IF(DMP)  289 r 289 r 286  B1520 

286  WRITE (6»287)  AO  (Jr I) r BB ( Jr  I ) r CB  < J» I) rDB ( Jr  I ) rJrlrA(L) rB(L) rC  (L) r  D(  B1530 

!L)fL  B1540 

287  FORMAT (1H0*8HAB (Jr  I ) = r 1PE12 . 5 r 2X r 8HBB < J r I ) = r 1PE12. 5 r 2X r 8HCB ( J r I ) = r  B1550 

HPE12.5r2Xr6HDB(JrI)=,lPE12.5r2Xr2HJ=rI3,2Xr2HI=rI3/lXr5HA(L)=rlPE  B1560 
212.5 r2Xr5H8(L)=r 1PE12 . 5  r  2X  r  5HC  (L )  =  r  1PE 12 . 5r  2X  r  5HD  ( L)  =  r  1PE12 . 5  r  2X  r  2  B1570 

3HL=rI3)  B1580 

289  L=L+1  B1590 

290  CONTINUE  B1600 

300  CONTINUE  B1610 

NPFT=L-1  B1620 

RETURN  B1630 

END  B1640 
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* IBFTC  PRP  COOOO 

c  this  subroutine  dftfrminfs  the  physical  properties  of  thf  room 

C  HEAT  SHIELD  STRUCTURE  COORD 

SUBROUTINE  PROP  C0030 

C  C0040 

DIMENSION  TlTl E ( 1? ) ,HFADNG(12) , XIDNT (12.12) ,TKC(20) ,  XKC(20> ,  COOSO 

1CPC ( 20 ) ,  TKV ( 20 ) , XK V ( 20 ) »  TOPV ( 20 ) . CPV ( 20 ) ,TIMF(30o) .OCON(300) ,  00060 

2ORA0  (  300  )  , VEL  <  300 )  ,Xf\lPM(  )  2)  ,NKPR(  12)  ,NCPR( 12)  ,TXK  (20, 12.)  ,XK(20» 12)  C0070 

3. TCP ( 20 , 12) ,CPX ( 20, 1 2) ,PH0BX(12) ,XBM() 2) ,FMFB( 12) »EMRR( 1?) ,HXX( IP)  COO AO 

4  » G  APX ( 12 ) ,FTfST(12> ,flTEST(12) .TEMPI (200) ,TX1  ( 20  0 ) , TX2 ( 200 ) .  C0090 

5TX2T (10*12) .TltLl (200) ,TUL2(200) ,HX ( 50 ) . Tw ( 50 > , IR ( SO ) , TR1 < SO ) ,  COlOO 

6TR2 ( 50 )  *  TUL  ( 50 )  ,  IFM(SO)  , TY ( 200 ) , A ( 200 >  ,R(200)  ,C  (  200  )  , D ( 2n0  )  ,  COHO 

7R ( SO ) ,RH0(50) ,CP(50)»DXP(12)»XKB(10,12)»CPB (10.12) *XMnG(SO>  *  CO  120 

flYK (50) *  AH (10. 12) *BB( 1 0, 12) *CR(10.12) . DB ( 10.12) .SR (10, 12) .  C0 130 

9RPK10. 12)  *RR2(10.12)  ,H(  12)  .S(50)  .NPM(12)  C0140 

DIMENSION  TTl)L(SO) .RH0Y1 (SO) » RH0Y2 (50 ) .DRHO(SO) .TCPC(PO)  C0150 

C  CO  iso 

COMMON  TKC,XkO»TCPC,CPC.TKV.XKV.TCPV.CPV,XNPM,RHOBX,XPM.FMBR,  CO 170 

1fMFB,NKPB.NCPR.TXK,Xk,TCP,CPX,NPM.GAPX.FTFST.BTFST.TEMDI,TX1 ,  C01A0 

2TX2,TX2T.TUL.TUl 1.TUL2.IR.IR1 ,IP2.A,B,C.D.S.R,AP,RB.CR.DR.SP.  00190 

3RP1.RP2. TY.RhOYI .RH0Y2.XMDG, RHO, CP. YK,XKr,CPB,DXR.DT,XLOsT,  00200 

4TABL.TCHAP.TRFC.RH0V , RHOC .FBLOW • fMV , EMC .H300 . NKO , NCPC . NKv »  NCPV .  0021  0 

5NP,NMR,NPRS,NPF.TEST2,  TEMPI,  TX0,TFNV,HEN\/,FFNV,0|..OSS,TLlM.  TINT  C0220 

COMMON  I1.I2.T3.I4.IR. 16 ,OIN , TNT . DX , XMT , TL . VL . BL, DMP, FRP 1 ,EPR2.  00230 

lFRR3.ERR4.HV, VPT.CHARK .CHARC, ABLK. ARLC.XmDC.H  00240 

C  C02S0 

KTNT=TNT  00260 

DO  170  I=KINT,NP  00270 

10  IF ( IR ( I >  >  12,12.100  C02A0 

12  TUL < I)=AMAX1 (TX1 ( I ) ,TX2< t ) )  C0290 

TF(TUL( I )  .LF.TAPL)  GO  TO  20  00300 

TP(I)=1  00310 

fiO  TO  lOu  00320 

20  TF ( I 1—1  25.2  . 2l  00330 

21  TF(Il-NKV)  ,-'2.:>2.25  00340 

22  IF (TX2(I )  —  TK  v  < II ) )  35,55,30  00350 

25  WRITE ( 6 . 25)  TX2(I)  00360 

26  FORMAT (1H0.A7H  THF  RANGE  OF  ONE  OF  THF  ARLATION  PROPERTY  CURVF  FTT  00370 

IS  WAS  EXCEEDED  AT  A  TFMPFPATURE  OF  .1PE12.5)  C03A0 

FRR2=1.0  00390 

GO  TO  355  00400’ 

30  I 1 si 1+1  00410 

GO  TO  21  00420 

35  TF(TX2(I) — TKV ( I 1— 1 ) )  40.SS.50  00430 

40  11=11-1  00440 

GO  TO  20  C04RO 

SO  YK ( I ) =XK V ( 1 1  —  1 )  +  ( <XKV( 11 )-XKV(Il-l))/[TKV(Il )-TKV(Il-1  ) ) )  00460 

1  * ( TX2 ( I ) — TKV ( 1 1— 1 ) )  00470 

GO  TO  60  C04R0 

55  YK ( I ) =XK V (II)  00490 

60  IF ( 12-1 )  25.2S.6l  00500 

61  TF(I2— NCPV)  62.62,25  00510 

62  IF (TX2  < I ) -TCPV ( 12 ) )  70, PS, 65  00520 

65  12=12+1  00530 

GO  TO  61  00540 

70  IF(TX2(I )-TCPV( 12-1) )  75.P5.flO  005SO 

75  12=12-1  00560 
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GO  TO  60 

C0570 

AO 

CP ( I ) =CPV  < 12-1 )  +  ( <CPV( I2)-CPV( 12-1))/ (TCPV(X2)-TcPV( 12-1) ) ) 

C05B0 

1* ( TX2 ( I ) -TCPV ( 12-1 ) ) 

C0590 

GO  TO  90 

C0600 

A5 

CP(I)=CPV(I2) 

C0610 

90 

RHO(I)=RHOV 

C0620 

GO  TO  170 

C0630 

100 

TUL ( I ) =AMAX1 (TUL  < I ) »TX2 (X) ) 

C0640 

XF(TUL(I)-TChAR)  110,110,115 

C0650 

110 

RHO ( I ) =RHOV+ { RHOV-RHOC ) * ( <  TUL  < I ) -T ABL ) / <  T ABL-TCH AR ) ) 

C0660 

YK ( I ) =CHARK+ ( ABLK-CHARK )  *  (  < RHO ( I ) -RHOC ) / ( RHOV-RHoC ) ) 

C0670 

CP ( I ) =CHARC+ ( ABLC-CHARC ) *  < ( RHO ( 1 1 -RHOC ) / ( RHOV-RHoC ) ) 

C06A0 

GO  TO  170 

C0690 

115 

XF(VPT)  116,116.117 

C0700 

116 

TTUL(I)=TUL(X) 

C0710 

GO  TO  120 

C0720 

117 

TTuL ( X )  =  TX2 ( I ) 

C0730 

120 

XF ( 13-1 )  25,26,121 

C0740 

121 

XF ( I3-NKC )  122,122,25 

C0750 

122 

XF(TTML(1)-TkC(I3) )  124, 135,123 

C0760 

123 

13=13+1 

C0770 

GO  TO  121 

C07A0 

124 

TF(TTUL(I)-TKC<I3-D)  125.135.130 

00790 

125 

T3=I3-1 

COAOO 

GO  TO  120 

COftlO 

130 

YK ( I ) =XKC ( 13-1 )+( (XKC( 13) -XKC( 13-1))/ (TKC(I3)-TKC( 13-1 ) ) ) 

00820 

1*(TTUL(I)-TKC(I3-1)) 

C0B30 

GO  TO  140 

C0A40 

135 

YK ( I ) =XKC ( 13 ) 

COB50 

140 

XF ( 14-1 )  25,26,141 

C0A60 

141 

TF ( I4-NCPC  )  142,142,25 

COA70 

142 

XF  ( TTI IL  ( I  )-TcPC  ( 14 )  )  150,165.145 

00880 

145 

14=14+1 

C0A90 

GO  TO  141 

C0900 

150 

)F<TTULU)-TCPC(I4-1  )  )  155. 165.160 

00910 

155 

14=14-1 

C09?0 

GO  TO  140 

00930 

160 

CP ( I ) =CPC ( 14—1 )+( (CPC (I4)-CPC(X4_1) )/{TCPC(X4)-TcPC (14-1 ) ) ) 

00940 

1* ( TTUL ( I ) -TCPC ( 14-1 ) ) 

00950 

GO  TO  166- 

00960 

165 

CP(I)=CPC(X4) 

00970 

166 

RHO(I)=RHOC 

009B0 

170 

continue 

00990 

01000 

oftermination  of  propfr  back-up  shield  material  property 

01010 

01020 

01030 

01040 

rin5n 

C1060 

01070 

oioao 

01090 

cuon 

rmr> 

01120 

.I2.71H  BACKUP  5TPllC  01130 


no  3oo  i=i.nmr 

nXB(I)=X6M(I)/(  (XNPM(I)-I  .0*12.0) 

LKP=NKPB(I) 

LCP=NCPB(I) 

NN=NPM(I) 
no  280  J=1»Nn 

200  IF ( X5-1 )  203.203.201 

201  IF(I5-LKP)  202. 202.203 

202  IF(TX2T<J,I)-TXK(I5,I) )  206. 220. 205 

203  wPITE (6.204)  X.TX2T(j,I) 

204  FORMAT (1H0.32H  THF  RANGE  OF  ONE  OF  THF  NiiMBFR 
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1TURF  PROPERTY  CURVE  FIT*;  WAS  FXOFFDFD  AT  A  TFMPFRATURF  OF  . 1PF12.5  01140 

2)  Cl  1 50 

FRR2=1.0  01160 

00  TO  35b  01170 

205  15=15+1  Cl 1 80 

00  TO  201  01190 

206  IF(TX2T(J»I )— TXK ( T5- trill  210, 220.215  C1200 

210  T 5=15-1  C1210 

00  TO  200  C 1220 

215  XKH<J,I)=XMIS-l,I)  +  <  (  XK  (  T5.  T  )— XK  (  15—1 , 1  )  )  /  < TXK  (  T5.  T  ) -TXK  <  15-1  .11  >  C1230 

I) *(TX2T<J,I)-TXK(I5-1,I) )  ClPttO 

00  TO  230  C1250 

220  XKB(J,I)=XK<IS,I)  C1260 

230  TF ( 16-1 )  203,203,231  C1270 

231  xF ( I6-LCP)  232.232. 2n3  C1280 

232  TF(TX2T(J.I)-TCP(T6.T) )  234,205.233  C1290 

233  15=16+1  C1300 

00  TO  231  C1310 

234  TF(TX2T(J,I)-TCP(I6-1 ,1) )  235,24s, 240  C1320 

235  16=16-1  C1330 

00  TO  230  C1S40 

240  CPB<vJ,I)=CPX< 16-l.I)  +  <  <CPV<l6,I)-CPX<T6-1  .1)  )/(TcP(l6,I)-TCP(T6-1  .  C1350 

II)  ))*(TX2T(J,T)-TCP<I6-1,T))  C1360 

00  TO  280  C1370 

245  CPB(J,X)=CPX(16,I)  C1380 

280  CONTINUF  01390 

15=2  C1400 

16=2  C1410 

300  CONTINUF  C1420 

310  IF(OMP)  355,365. 320  C1430 

320  WRITE  <  6, 330 )  01440 

330  FORMAT ( /IX. 32H  PROPFRTIFS  OF  ABLATION  MATFRIAL/I  C1450 

WRITE ( 6 »  335 )  C1460 

335  FORMAT ( /5X , 5HYK ( I ) » 9x »  5HCP ( I ) »9X • 6HRH0 1 1 ) / )  C1470 

wP ITE ( 6 , 340 )  <YK(I) ,cP(I) ,RHO ( I ) , 1=1 , NP )  C1480 

340  fORMAT(2X.1Pfi2.5,2X,1PF12.5,2X»i PE12 . 5 )  01490 

WRITE ( 6 »  345 )  C1500 

345  FORMAT ( // 1 X , 32H  PROPFRTIFS  OF  BACK-UP  STRUCTURE/)  C1510 

WRITE ( 6 , 347 )  C1520 

347  FORMAT (/5X.8hvKR( J. I) , 7X , 8HCPF < J , I ) ,7Y,8mRH0BX(I) , 7X, 7HFmFB ( I ) , 8x ,  C1530 

17HFMBR(I) ,9X,6HDXB(I)/)  C1540 

no  350  1=1, Nmp  01550 

KLsNPMII)  C1560 

no  349  J=1 ,KL  C1570 

WRITE (6, 348)  VKB(J.I) ,CPR(J»D .RMOBX <T) ,FMFR ( I > , fMBR (T> , OXB ( T )  C15S0 

348  F0RMAT(3X,1Pf12.5,3X. 1 PF1 2. 5 , 3X . 1 PE12 . 5 , 3X , 1 PF12 . 5 , 3X , 1PF1 2 . 5 , 3X. 1  C1590 

1PF1P . 5 )  C1600 

349  CONTINUE  01610 

350  CONTINUE  C1620 

355  rFTURM  01630 

FND  C1640 


sibftc  ABL  nooon 

c  this  subroutine  determines  thf  mass  flow  patf  from  thf  nooio 

c  ablating  NonES  no020 

SUBROUTINE  ABLATE  D0030 

c  nooun 

DIMENSION  TITLE! 12) »HEA0NG(12) .XtONT( 12.12) .TKCI20) .XKC(20> ,  00050 

lrPC<20).TKV<20) .XKV(20) » TCPV (20),CPV(20)»TIME(300) .QCON(300) ,  00060 

2oPAD(300) »VEL(300) »XnPM(12)  »  NKPP  ( 12 ) ,  NCPB !  12 )  .  TXk  « 20 , 1  2 )  ,  XK  (  20 , J  2 )  00070 

3.  TCP (20 .12) .CPX< 20.12) .RH0BXO2) ,XBM(12> ,FMFB(12) .EMBRI12) .HXX(l2)  00080 

4, GAPX(12> .FTesT(12) .BTEST<12) ,TEmOI<200) ,TX1 (200) .TX2<200)  .  00090 

5TX2T(10»12)»TUL 1(200) ,  TUl.2  ( 200 )  .  HX  ( 50 ) , Tw ( 50 ) , IR ( 50 ) , TR1 ( SO ) ,  001 00 

6lP2(bO) » TUL (SO ) • I EM (SO) .TY(200) . A (200) »B(200),Ct200).n(200),  00110 

7r ( 50 ) » RHO (50 ) , CP ( 50 ) ,0XP(12) ,XKB<10.12> »CPB(10*1?) . XMOG(SO) .  D0120 

8YK (50) .Afclt 10,1 2) ,RB< 10.12) »CR{ 10, 12).0B( 10.12 ) .SR( 10,12) ,  001  SO 

PRPK10.12)  *  Rr2  (10.12)  »  H  ( 1  2  )  .5 < 50 >  » NPM < 12)  00140 

DIMENSION  TTuL ( 50 ) *RhOY1 (SO) ,PHOy2(50) »0RH0(50) .TCPC(20)  00150 

C  00160 

common  TKC.XkC.TcPC.CPC.TKV.XKV.  TCPV.CPV, XNPM.RHoBX.XRM.  EmBB  >  00)70 

1fMFB»NKPB.NCPB>TXK»Xk»TCP»CPX,NPm*GAPX»FTEST»RTEsT»TEMDT,TX1 ,  DO  180 

2TX2f TX2T» TUL, TUL 1» TUL2.TR *IR1,IR2. A, B,c»n, S#R,AB,BB »CB»DB* SB.  00190 

3rP1,RB2#TY.RHOY1,RH0y2*  XMPG ,RHO »CP.YK»XKR,CPB,DXR,DT, VLOST ,  00200 

4TABL,TCHAR.TRFC»RHOV,RHOC.FRLOW.EMV.EMC»H30n,NKC.NCPC.NKV.NCPV»  00210 

5NP,NMR»NPFS»NPF.TFST3,TFMPI .TXO.TENV.HENv, fenv, GLOSS, TLIM. TINT  00220 

COMMON  1 1 , 12, T3, 14, IS, 16,0 IN, INT , OX , XMT . TL . VL , BL ,OMP , ERRl , ERR2 »  00230 

lERR3,ERR4,HV,yPT,CHARK .CHARC. ABLk. ABLC.XMOC.H  00240 

C  00250 

xMT=0.0  00260 

LTNTfTNT  00270 

kI=NP  00280 

IF(DMP)  ri.8.3  00290 

3  WRITE (6 . 5)  00300 

5  fORMAT(//1X.2PHMASS  FLOW  FROM  ABLATING  NODES//)  00310 

8  DO  200  KKT=LlNT,NP  00320 

TF(IRKKI))  11.11,12  00330 

11  TF(TX1 (KI) .LE.TABL)  SO  TO  9  00340 

12  TUL1 (KDaAMAX) (TULI(KT) .TYl(KT) )  00350 

TP1(KT)=1  00360 

GO  TO  20  00370 

9  IF(TX1 (KI)-TaPL)  10,10.20  00380 

10  RHOYKKI  )=RHow  00390 

GO  TO  50  00400 

20  TF(TUI.KKI)-TCHAR)  40,30,30  00410 

30  RHOYKKI  )=RH0C  00420 

GO  TO  50  00430 

uo  RHOYl ( KT ) =RH0V+ (RHOV-RHOC ) * ( (TULl (KI)-TARt  )/(TAP|  -TCHAR) )  00440 

50  TF(IR2(KI)>  52.52.54  n0450 

52  TE(TX2(KI). LE.TABL)  GO  TO  56  00460 

54  TUL21KI )=AMAX1 (TUL2(KT) .TX2(KT) )  00470 

IP2(KT)=1  00480 

GO  TO  70  00490 

56  TF(TX2(KI)-TAPL)  60.ftO.7O  00500 

60  RH0Y21KI )=RHOV  00510 

GO  TO  95  00520 

70  TF<TUL2(KT)-TCHAR)  9o, 80,80  00530 

80  RH0Y2(KI )=RHOC  00540 

GO  TO  95  00550 

90  RH0Y2 (KI ) =RHOV+ ( RHOV-RHOC ) * ( (TUL2 ( KI >-T ARL > / ( TABL-TCHAR ) )  00560 


95  ORHO(KI)  =  (  <RHOYl<KI)..RHOY?<Kn  )/OT)*DX 
TF(KI-NP)  97,96>96 

96  nRH0(KI)snRHn(KI)/2.n 
GO  TO  98 

97  tF(KI-INT)  96,96,98 

98  IF < DRHO ( KI ) )  110.120,120 
110  nRHO(KI)=0.0 

120  XMT=XMT+DRH0(KI) 

XMDG<KI )SVMT 
TF(DMP)  100.190,150 

150  wPITE (6 » 160 )  XMDG (KI ) , DRHO < K I ) , RhOY2 ( KI ) , RHOY1 ( K I ) 

160  F0RMAT(lX,5HXMDG=,lPpl2.5,2X,5HDRH0=,lPFl2.5.2X,6HRH0Y2=,lPEl?.5.2 
lX, 6HRH0Y1S, 1PF12 . 5 ) 

190  kI=KI-1 
200  CONTINUE 
RFTURN 
FND 
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SIBFTC  0X10 

THIS  SUBROUTINE  CALClll.ATFS  THF  HfATING  RATE  DUE  TO  COMBUSTION 
IT  IS  ASSUMED  THAT  OXYGEN  AND  CARBON  REACT  TO  FORM  CO  ONLY. 

SUBROUTINE  OXTDAT(XMDO,OOXID) 

C 

OOX  ID=XMDO*4000. 0/36(10.0 

(30X10=0,0 

RETURN 

END 
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SIBFTC  SWUFT 

C  THIS  SUBROUTINE  DFTErMINfS  THF  FORWARD  TIME  STEP  TEMPFRATURFS 

C  RY  SOLVING  The  tri-diagonal  matrix 

SUBROUT I NF  SWUFT { A , R , r  * n . T » N , nMP ) 

dimension  A<?no),R<2nn>*c<200).D(2oo>»T(?oo),rp<200),np<?oo) 

CP(1»=C<H/B(1 > 
nPCl)=D(l)/B(1) 
no  ion  i=?»n 

CP(I)=C(I)/(BCI)-A(I)*CP(I-1)> 

DP< I )=(D( I )-A( I)*DP( T-l) )/<B<  T  >-ACI)*fP( T-l) ) 

100  CONTINUE 
T(N)=DP(N) 

NM1-N—1 

no  2on  j=i »nmi 

T=N-J 

T(I)=nP(X)-CP(I)*T(I+1) 

200  CONTINUF 

TF(DMP)  300  f  300 » 250 
2S0  WRITE (6*260 ) 

260  FORMAT(//lX »43HC0FFFi CIFNTS  CALCULATED  By  SURROUTINF  SWUFT//) 
wR ITE ( 6  *270 ) 

270  FORMAT! 6X. 5HCPII) »10X.5HnP(I) ,10X»4HT(I>/) 

WRITE (6 * 275)  <CP< I ) ,nP< I ) *T < 1 ). 1=1 *N) 

275  fORMAT!PX*1PE12.5*2X, 1 PFl 2.5 *PX *  1 PE12 . 5) 

300  rFTURN 
FND 
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sibftc  rec  ooooo 

oooio 

this  subroutine  determines  the  front  facf  location  and  char  mass  goopo 

REMOVAL  RATE  G0030 

G0040 

SUBROUTINE  RFCESS ( XMnC *  XI  OST ,TRFC . DT , RHOC . TS . SR .TX2 . NPEC , NRS ,FRRs .  00050 

lsXO  .SOOT  .  OMP)  G0060 

C  GO 070 

DIMENSION  TS ( 50 ) »  SR ( 50 )  G0080 

TF ( TX2-TRFC )  10. 20.20  G0090 

10  XMDC=O,0  G0100 

XLOST =0  ,  Cl  GO  110 

SDOTrO.O  G0120 

GO  TO  60  00130 

20  IF ( NRS-1 )25*25*21  00140 

21  IF ( NRS-NRFC  )  22.22*25  G0150 

22  TF ( TX2-TS ( NRS ) )  32.4n,30  G0160 

25  wPlTE (6  *  26 )  TX2  G0170 

26  FORMAT ( 1H0  *  75H  ThE  RANGE  OF  THE  SURFACE  RECESSION  TA0l  E  WAS  EXCFFD  00180 

1FD  AT  A  TEMPERATURE  OF  .1PE12.5)  G0190 

FPR5=1.0  G0200 

GO  TO  60  0021 0 

30  NPS=NRS+1  00220 

GO  TO  21  O0230 

32  TF(TX2-TS(NRS-1) )  34,40.36  O0240 

34  NPS=NRS-1  00250 

GO  TO  20  00260 

36  SX=SR(NRS-1)+((SR(NRS)-SR(NRS-I)l/<TS<NRS>-TS(NRS-I)n  00270 

1* ( TX2-TS ( NRS-1 ) )  00200 

GO  TO  50  G0290 

40  SXrSR(NRS)  G0300 

SO  VI OST=300.0*SX*DT  00310 

XMDC=(XLOST*RHOC)/DT  00320 

SDOT=SX*300,0  00330 

IF(DMP)  60.60.52  00340 

52  WRITE ( 6. 54 )  SX. XlOST , YMDC  G035O 

64  FORMAT  1 1H0 . 3HSXS . 1PE 12.5. 3X  *6HXL0ST=* 1PE 12,5, 3X» 5HXMDC- . 1  PEI  2,5)  00360 

60  RFTURN  G0370 

END  00300 
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SI8FTC  TEMP  HOOOO 

THIS  SUBROUTINE  DFTERMINFS  THF  INITIAL  TFMPFRATURF  DISTRIBUTION  HOOIO 

IN  THF  HEAT  SHIFLO  STRUCTURE  H0020 

SUBROUTINE  TEMPD  M0030 

C  H0040 

DIMENSION  TlTl  E  ( 1 ? ) .HEAONG ( 12 ) »XlDNT (12,12), TKC(20),XKC(20) »  HOOSO 

1CPC (20) » TK V ( 20 ) » XKV ( 20 ) #  TCPV  ( 20 ) »CPV ( 20 ) .TIME (300) »GCON(300) .  H0060 

20PAD ( 300 ) »  VEL (300)»XnPM(12) ,NKPR( 12) »NCPB(12) , TX« ( 20 . 1 2 ) , XK ( 20 » 12 )  H0070 

3.  TCP (20 ,12) ,CPX(20»l2) # RH08X ( 12 ) . XBM< 1 2 ) , FMFB ( 12 ) ,EMBR ( 1 2) »HXX ( 12 )  HOOAO 

4 .  GAPX  (12)»FTEST(12)»RTEST(12)»TFmDI(200)»TX1(200)»TX2(200) .  H0090 

5TX2T(10»12)»TUL 1(200), TUI.  2 ( 200 ) » HX ( 50 ) »Tw(50) *IR(50) .IR1 (50) ,  HOI  00 

6rP2(50) *TUL(SO) . IEM (SO) , TV (200) . fl ( 200 ) »B(200) *C(200) *0(200) ,  HOI 10 

7r ( 50 ) .RHO ( 50 ) ,CP(50).OXB(12).XKB(10.12) .CPB ( 10 . 12) , XMDG( SO ) ,  H0120 

8YK(50) »AB( 10,1 2), BB( 10,12 )»CR( 10.12 ) ,DB( 10 ,12) .SR (10,12 ) ,  H0130 

9RBK10.12)  ,RB2(10.12)  ,H(12)  .S (50  )  ,NPM{  12 )  H0140 

DIMENSION  TTlJl.(50)  .RHOYl  (SO)  ,RHOy2(50)  ,DRHO  ( 50 )  ,TCPC  ( 20 )  H0150 

C  HO 160 

COMMON  TKC.XkC.TcPC. CPC, TKV.XKV, TCPV, CPV.XNPM.RHoBX.XRM.FMBB,  H0170 

1FMFB,NKPB»NCPR,TXK,XK,TCP.CPX,NPM,GAPX,FTFST,RTEST»TEMDI.TX1,  HOI AO 

2tX2,TX2T,TUL.TUL1»TUl2, IR , IR1 , IR2, A , B, C »D» S.R , AB.BB .CB.DB. SB ,  H0190 

3RB1.RR2.TY. RHOYl, RH0Y2.XMnG, RHO, CP, YK.XKR, CPB, DXb.DT.XLOST.  H0200 

4TABL , TCHAP , TRFC ,RHOV , RHOC  »FBLOW,fMV , EMC  »H300 » NKC. NCPC » NKv ,NCPV,  H0210 

5NP , NMR , NPBS , NPF » TFSTp, TEMPI »TXO » TFNV , HEN V , FFNV , OLOSS , TL IM , T I NT  H0220 

COMMON  11,12.13,14,15, 16 , ©IN, TNT ,DX, XMT »TL » VL , BL ,OMP, FRPl , ERR2 ,  H0230 

1FRR3.FRR4.HV,VPT,CHARK.CHARC. ABLK.ABLC.XMDC.H  H0240 

C  H0250 

*=0.0  H0260 

IF ( TEST2 )  300, 100. 200  H0270 

100  DO  150  L=1»NPF  H02AO 

TXl (L)=T£MPI  H0290 

TX2(L)=TEMPI  H0300 

TUL1(L)=TXKL)  H0310 

TULP(L)=TX2(l)  H0320 

TFMDI(L)=TEMPI  H0330 

150  CONTINUE  H034O 

DO  160  1=1, NmR  H0350 

jNsNPM ( I )  H0S60 

DO  155  M=1.JN  H0370 

TX2T(M,I)=TEMPI  H03A0 

155  CONTINUE  H0390 

160  CONTINUE  H0400 

GO  TO  320  H0410 

200  no  220  L=1.NP  H0420 

TEMDI (L)=TXO+( (TENV-tXO) /TL)**X  H0430 

TXl (L)=TEMDI (L)  H0440 

TX2(L)=TX1 (L)  H0450 

TUL1 (L)=TX1 (L)  H0460 

TUL2 (L ) =TX1 (L )  H0470 

XSX+DX  H04A0 

220  CONTINUE  H0490 

L=NP+1  H0500 

nO  270  1=1 »NMR  H0510 

KJsNPM(I)  H0520 

DO  250  Jsl.Kj  H0530 

TEMDI (L)=TX0+( (TENV-TXO) /TL)**X  H0540 

TX1(L)=TEMDI(L)  H0550 

TX2 (L) sTEMDI ( L )  H0560 
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TX2T(J,I)=TEMni<D  H0E7O 

X=X+DXB(I)  HOSao 

L=L+1  H0S90 

2E0  CONTINUF  H0600 

X=X+(GAPX<I)/12.0)  H0610 

270  CONTINUF  HOf.20 

GO  TO  3211  H0630 

AM  ARRRTTARY  TEMPFRATllRF  nlSTPIPllTION  CAN  BF  PEAn  IN  FROM  INPUT  HOfittO 

DATA  IF  TEST?  IS  A  NFGATTVF  NUMrFR  M0G50 

300  WRITE ( 6, 31 0 )  H06B0 

310  FORMAT(1HO*79H  THE  VALUF  OF  TFST ?  WAS  NFGATTVF,  SUBPOUTINF  TFMPO  S  M0fi70 

lHOllLD  NOT  HAVF  BEEN  CALLED.)  HOEflO 

ERR1=1.0  M0690 

320  pFTURN  H0700 

FND  H0710 


► 
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IBFTC  D0N2  TOOOO 

THIS  SUBROUTINE  DETERMINES  THF  TfMPFRATUrF  OF  POINTS  A  FIXED  T0010 

DISTANCE  FROM  A  RFFERFNCF  PLANE  FROM  THE  TEMPFRATIJRES  CALCULATED  T 005>0 

in  a  varying  thickness  iooso 

10040 

SUBROUTINE  DoM2  (  XLOST  »XARRAY«TARRAY*NAt  XnODE  >  TEMP  »XNOOEVfKK»  XLSTv  >  lOOSO 
lnx)  10060 

C  T0070 

DIMENSION  XARRAY(SO) ,TAPRAY(50) *VNODE(50>  , TEMPI  So > *XNODFv(50)  iooro 

C  10090 

K=0  T0100 

nXT=0.0  T0110 

no  loo  1=1 »na  10120 

IF < XLSTV, LE.DXT)  GO  TO  ISO  T 0 1 30 

K=K  +  1  T  0 140 

100  nXT=DXT+DX  T0150 

ISO  kK=NA-K  T0160 

VK=K  TO  1 70 

XNODEV< 1 )=XLSTV  TOIflO 

TFMP ( 1  ) =T  ARRAY ( 1 )  T019O 

no  200  1=1 rKK  T0200 

yNODE ( I ) =XK*nx-XLSTV  T0210 

call  DISCT3 ( XNODE ( I  )  ,  XARR  AY  *TARRAY»NA» TFmP  1 1  +  1 )  )  T0220 

XNODEVI 1+1 )=XK*OX  T0230 

200  VK=XK+1.0  T0240 

RETURN  T02SO 

FND  T0260 
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SIBFTC  UINTRP 

joono 

FUBROUTInF 

UINTRPTX.XTBI ,Y*YTPL#N* J) 

J001O 

niMFNRION 

XTHI.  (50)  >YtRL<50) 

J0020 

T=J 

J0030 

TFCI.fiT.N. 

0R.T.LT.2)  T  =  2 

J0040 

10 

TF(XT«U( I- 

l).TE.X.ANn.X.lF.XTPL(T))  GO  TO  40 

J0050 

TF  <  X .AT . XTBL ( T ) 1  GO  TO  30 

J0060 

20 

T  =  I-1 

J0070 

TFII.AE.21 

GO  TO  10 

JOOAO 

T  — 2 

J009O 

AO  TO  40 

J010O 

SO 

T  =  I  +  1 

J0110 

TF(I.LE.N) 

Go  TO  10 

J0120 

t=n 

J0130 

UO 

FRACT=(X-XTBLd-l)  )/(XTPLd)-XTPL(I-l)  1 

JOl  40 

Y=YTBL  <1-1  )  +  <YTR|_<Il-YTPl  <1-11  )*FRACT 

JO  150 

rftukn 

JO  1ft  0 

fnd 

JOl  70 
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SIBFTC  TSOT 

FUBROIITINF  ISOTHM (DEPTH. TEMP, PONn.N.ANS) 

DIMENSION  DEPTH! 1)»TPMP(1  ) 

AWS=-1 . 
k=N-1 

no  ion  t=i.k 

tf ( temp ( i ) -Bond )? » l »t 

1  ANS=DFPTH(I) 

00  TO  100 

2  lF(TEMP(i+l)_ROND)10n,lnn»4 

4  AMS=DFPTH ( I  + 1 1 - ( TEMP ( I  + 1 ) -BONO ) * ( DEPTH ( 1  +  1  > -DEPTH ( I ) ) / ( TFMP ( I  +  H - 
lTFMP(T) ) 

00  TO  100 

3  TF(TEMP(I+1) -POND) s.ino.ino 

5  ANS=  ( TEMP ( 1 )-POND )*(DFPTH(I+1 )— DEPTH (!))/( TFMP(I )-TEMP( 1+1 )  ) +DEPTH 
1(1) 

lno  CONTINUE 

TF ( BOND. £0. TEMP (N) ) ANSzDFPTH (N ) 

RETURN 

FND 
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sibftc  save  loooo 

FUBROIlTlNF  SAVE (SAVE l , SAVF2 »SAVF3»USE,NX1 , VALUE »DT . TFINAL » TIMF ,  L0010 

1THING)  l  0020 

OIMENSION  SAvFl ( 1 >  *  SAVE?  (t)#SAVF3(l)  1.0030 

USE=0.0  10040 

9AVEl(NXl)=VAt.UE  L0050 

NX2=NXl-l  L0060 

TF(NX2.FQ.0)NX2=3  1.0070 

SAVF2 (NX2 )= VALUE  LOOOO 

MX3=NX2-1  l 0090 

TF(NX3.EQ.0)nX3=3  1.0100 

SAVF3(NX3)=VA(  UF  1.0110 

TF((TTME.LT,( ?,*0T ) ) .OR. (TIMF ,GF. <TF INAL-3 »*DT) ) ) 60  TO  4  L0120 

GO  TO  (l»?*3),NXl  L  01 30 

1  TF( ( (ABS(SAVF?(1 )-SAvF2(?)) ) .LE. . 001 ) .OR . ( ARS ( SAvE2 ( 2 ) -SAVE2 ( 3) )  L0140 

1. IE. .001)160  TO  5  L0150 

TF( < (SAVE251) ,LT.SAVf?(?) ).ANO. ( 9A VF2 ( 2 ) .GT.SAVF? ( 3) ) 1 .OR. ( (SAVF?(  L01G0 
111.  T ,SA VF2 ( 2 1 ) .AND. (9AVF2(2) ,LT.SAVE2(3) 1 1 1 USE=sAVFp ( 2 )  10170 

5  THING=SAVF2<2)  L0100 

GO  TO  4  L0190 

2  TF(  < ( ABS(SAVF3(l)-SAvF3(?) > 1 .1  E . . 00 1 1 . OR . ( ABS ( SAvF3 ( 2 1 -S AVE3 ( 3 ) )  L0200 

1. IE, .001)160  TO  6  L0210 

TF<  (  ( SAVE 3(1)  . LT .SAVF3  (  2  )  )  .  AND.  ( SAVF3  (  2 )  .GT .SAVF^ (  3)  > ) .OR. (  (SAVFM  1.0220 

11 ) . GT ,SA VF3 ( 2 ) ) , AND, { FAVF3 ( 2 ) ,LT .SAVE3 ( 3 ) ) ) )USE=SAVF3 ( 2 )  l  0230 

6  THJNG=SAVF3(?)  L 0240 

GO  TO  4  1.0250 

3  TF(  (  (  ABSISAVF1  (1  )-SAvFl<?) )  )  .1, E. . 001 )  . OR .  ( ABS ( SAvFl  ( 2 ) -S AVE1  ( 3)  )  L0260 

1. IE. .001)) GO  TO  6  L0270 

1F( ( (SAVE1 (1) .LT.SAVF1 (2) ) . AND. (9AVE1 ( 2) .GT.SAVFl (3) ) ) .OR. ( (SAVE) (  L02B0 

11  )  .GT.SAVFl  (?)  )  .AND.  (FAVF1  (2)  .LT.SAVE1  (3)  )  )  )  USE=«;AVF1  ( 2 )  1.0290 

7  THING=SAVF1(21  1.0300 

4  MX1=NX1+1  10310 

tF<NX1,F0.4)NV1=1  10320 

RFTURN  L0330 

FND  L0340 
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:  PISCT3 

Mooon 

SUBROUTINE  DisCT3(XA,TABx»TARY«Ny» ANS ) 

MOOIO 

DIMENSION  TABV(l) .TABY( J ) 

M002fl 

CALL  OISSFR (XA»TARX#l*NYt?*NN) 

M0030 

NMN=3 

M004d 

CALL  LAGkAN<XA»TARX<nN) ,T ABY < NN ) , NNN, ANS > 

MOOSO 

RFTURN 

M0060 

fnd 

M0070 

74 


SIBFTC  DISS  WOOOO 

SUBROUTINE  DISSER  (XA»TAB> I »NX» ID»NPX)  W0010 

DIMENSION  TAB <2000 )  W0020 

C  DIMENSION  TAB (2000 )  W0030 

NPT=ID+1  W0040 

NPB=NPT/2  W0050 

NPU=NPT-NPB  W0060 

IF  (NX-NPT)  10>5,10  W0070 

5  NPX=I  W0080 

RETURN  W0090 

10  NLOW=I+NPB  W0100 

NUPP=I+NX-<NPU+1)  W0110 

DO  15  I I=NLOW  *  NUPP  W0120 

NLOC=II  W0130 

IF  < TAB < 1 1 ) “XA >  15»20»20  W0140 

15  CONTINUE  W0150 

NPX=NUPP-NPB+1  W0160 

RETURN  W0170 

20  NL=NLOC-NPB  W0180 

NU=NL+ID  W0190 

DO  25  JJ=NUNU  W0200 

NDIS=JJ  W0210 

IF  (TAB<JJ)-TAB(JJ+1) )  25#30»25  W0220 

25  CONTINUE  W0230 

NPX=NL  W0240 

RETURN  W0250 

30  IF  (TAB(NDIS)-XA)  40»35»35  W0260 

35  NPX=NDIS-ID  W0270 

RETURN  W0280 

40  NPX=NDIS+1  W0290 

RETURN  W0300 

END  W0310 
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sibftc  lagr 

T0000 

SUBROUTINE  LAGRAN  ( XA , X > Y »N# ANS) 

T0010 

DIMENSION  X (200 ) *  Y (200 ) 

T0020 

c 

DIMENSION  X (200 ) * Y (200 ) 

T0030 

SUM=0.0 

T0040 

DO  3  1  =  1  *  N 

T0050 

PROD=Y(I) 

T0060 

DO  2  J=1*N 

T0070 

A=X(I)-X(J) 

T0080 

IF  (A)  1»2»1 

T0090 

1 

B=(XA-X(J) )/A 

T0100 

PROO=PROD*B 

T0110 

2 

CONTINUE 

T0120 

3 

SUM=SUM+PROD 

T0130 

ANS=SUM 

T0140 

RETURN 

T0150 

END 

T0160 

SIBFTC  MORE  N0000 

DIMENSION  TITLE (12)  »X(2000)»Yl(2000)*Y2(2000)  ,Y3{2000>  .Y4(2000)  NOOlO 

REWIND  11  N0020 

READ(ll)  (TlTLE(I) .1=1.12)  N0030 

RFAD(11)X<1)  >Y1(1)  .Y2(l)  ,Y3(1) >Yu(l)  N0040 

y3(1)=Y3(1)*12.+Y1(1)  N0050 

Y4(1)=Y4(1)*1?.+Y1(1)  NOOAO 

1=2  N0070 

30  RFAD(11)X(I).Y1(I)»Y2(I),y3(I),Y4<I)  NOOftO 

TF(X(T)-5001. 510,20.20  N0090 

JO  Y3(I)=Y3(I)*12.+Y1(I)  NOlOO 

yMI)=Y4<I)*12.+Y1(I)  N0110 

T=I+1  N0120 

GO  TO  30  N0130 

20  NPL0T=I-1  N0140 

YM1=Y1(1>  N0150 

YM2=Y2<1>  N0160 

YM3=Y3 ( 1 )  N0170 

YM4=Y4(1)  N01R0 

DO  40  K  =  2  •  NPLOT  N0190 

TF  (Y1(K)  .GT.YM1)  YMi  =  yKKI  N0200 

TF  (Y?(K).GT.YM2>  YM?  =  Y2<K)  N0210 

IF  (Y3(K)  .GT.YM3)  YM.3  =  y3(K)  N0220 

IF  (Y4(K) .GT.YM4)  YMa  s  Y4(K)  N0230 

40  CONTINUE  N0240 

1000  FORMAT { 1H1 ,( 12A0 ) )  N0250 

CALL  ACCEND (X.Y1.Y2.Y3.  Y4 .NPLOT )  N0260 

XMAX=X<NPLOT)  N0270 

CALL  APLOT  (X.YI.XMaX, YMi .TITLE, NPLOT)  N02B0 

CALL  BPLOT  (X,Y2,XMAx, YM?, TITLE)  N0290 

CALL  CPLOT  (X.Y3.Y4.XMAX, YM3. YM4 .TITLE * Y1 )  N0300 

WRITE <6. 1000) (TITLE(T) »I=1  *12)  N0310 

WRITE(6»1001) (X(I>»Yl(I),Y2(I),Y3(I)»Y4<T),I=l * NpLOT)  N0320 

1001  FORMAT(5EPO.fl)  N0330 

WRITE (6# 1002 ) XMAX , YM t , YM2 , YM3 , YMu  *  NPLOT  N0340 

1002  FORMAT ( ///6H  XMAX=F1 0.4.SH  YM1=F10.4,SH  YM2=F10.a,5H  YM3=F10.4,5h  N0350 

lYM4=F10.4»2XftHNPLOT=I4)  N0360 

RFAD  (11)  (TITLE  (I), I  =  1.12)  N0370 

REA0(11)X(1).Y1(1) »Y2(1) ,Y3(1).Y4(1)  N0380 

1=2  N0390 

IF  ( X  ( l  )**5001 ,  )  30 .50  •  50  N0400 

50  WRITE (6. 1003) (TITLE< I) .1=1 .125  N0410 

1003  FORMAT (////12A6)  N0420 

RFTURM  N0430 

FND  N0440 
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SIBFTC  ACCEN  POOOO 

SUBROUTINE  AcCENO<X,YtA,B*C»N>  POOIO 

DIMENSION  X(1),Y(1).A<1).B(1).C(1)  PO02O 

K=1  P0030 

101  SMALL=X (K )  POOUO 

00  100  I=K»N  P0050 

DUMYsXd)  P0060 

SMALL=AMIN1(SMALL#0UMY)  P0070 

TF<SMALL.EQ.X(I>)INOeX=I  poobo 

100  CONTINUE  P0090 

X ( INDEX )=X (K )  P0100 

X<K ) -SMALL  P0110 

SAVEsY(K)  P0120 

Y(K)=Y<INDEX)  P0130 

Y< INDEX )=SAVE  P0140 

SAVEAsA(K)  P0150 

A(K )SA ( INDEX)  P0160 

A(INDFX>=SAVEA  P0170 

SAVEBrB(K)  POlflO 

B(K)=B(INOEX)  P0190 

R(INOFX)=SAVEB  P0200 

SAVFCrC ( K )  P0210 

C (K ) 2C ( INDEX )  P0220 

C( INDEX )=SAVEC  P0230 

K=K*1  P0240 

IF (K.FQ.N) RETURN  P02SO 

GO  TO  101  P0260 

FND  P0270 
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4IBFTC  /'(-'LOT  ooooo 

C.IIHROHTTMF  AP|  OT  (  X  •  Y ,  XI  TM ,  YLTM ,  Tl  Tt  E  ,  IP|  OT  )  OOOIO 

nIMFNRION  x< 30O) , yTTtlei i o) ,xtitlF( 10)  00020 

nIMFNRION  Tin  E( 1?) » Y(300) . ALONPYC7)  00030 

COMMON  /«PC  /  ALLOW ( 7) . Al ONGX (7 ) .NPLOT »7FR0, XMAX, TFIX  O004O 

DATA  CXTITLF(T)  ,I=1.10)/34H  TIME  (REC.)  /  OOORO 

OAT  A  (YTlTLF ( I > . 1=1. 10) /34H  RURFACE  RFCFRRIOm  (IN.)  /  00060 

7FKO=0.0  00070 

Al  LOW ( 1 ) =  R0.  00040 

At L0w(2)=1 00.  00090 

Al  LOW ( 3 ) =  260 ,  O0100 

Al  l.Ow(4)=ROO.  O01  1  n 

Al  L0w(5)=100o.  00120 

Al LOW(fe)=2bOO.  00130 

Al LOft(7)=EOOo.  00140 

MPL0T=IPLOT  O01R0 

no  lo  1=1.7  ooigo 

|I=I  00170 

IFIXLTM-  ALLOW ( 1 ) )  2 1 1 »  2  0  « 1  0  00140 

10  CONTINIJF  00190 

30  WRITE  (6,1000)  XLIM.YI  IM  00200 

1000  FORMAT ( ///77H  APL.OT  CANNOT  PF  DOmF  RECAHrF  FTTHPR  XLTM  FxrEFnFt)  RO  00210 
loo.  OR  YLTM  FvCFFOEn  R.  /AH  XI  IMrF 12 , R . R Y , 6H  YL T MrF 1 2 . ‘i  //  I9H  WF  O02?O 
2nOW  GO  TO  BPLOT  ///  )  00230 

PFTIIRM  00240 

20  vMAX=  ALLOW (IT)  O02RO 

1 F I X  =  T I  00260 

no  uo  1=1.4  00270 

Tlrl  00200 

1 F ( YL T M  *100.  -AlLOW(T)  )R0,G0.4n  00290 

UO  CONTINIJF  00300 

GO  TO  30  00310 

ro  ymax  =allow ( r t )  /ioo.  00320 

CALL  RSTFRM  O033O 

CALL  ORIOPN  (123. 1023, 24. 024, 1/3, ifl, 5, R)  00340 

CALL  PLOT  1  ( 1 , l.ZFRO, VMAx.ZFRO.YmAX.X, Y.wPLOT, 1 . 1M/)  O03R0 

Al  ONGX  ( 1  )=U.O  00360 

Al  ONGY ( 1  )=0,0  00370 

110  GO  1  =  1.6  00340 

CALL  LAPt.l  X  (  ALONPXd  )  .  1  )  00390 

CALL  LAREI  Y  ( ALONGY ( I ) , 1 )  00400 

Al ONGX ( T  +  1  ) =  ALONOX(I)  +.2*  XM  AX  00410 

GO  AI0NGY(I+1)=  ALONGY(I)  +.2*  YMAX  00420 

CALL  PRTNT(200.97R.l2.0.3«.XTTTLF )  00430 

CALL  PRINT(47,200.0. 1 2, 3H.YTITLF)  00440 

CALL  PHIinT(123. 1000. 12.0, 72, TTTI  F)  O04R0 

CALL  nMPbHF  00440 

RFTUKN  00470 

FNU  00440 
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■SIBFTC  PPL  OT  Roono 

FI iBPOl ITTnF  Ppi  OT  (X.y.XLTM.YL  tm.tTTLE)  pomn 

htmfnfion  x ( 3oo ) . y ( Xiin ) . ytitlf  u  n) . alongyi?) , ytttLE( in)  poopn 

I) T MFNF I OiM  TITIELI?)  P0030 

COMMON  /«PC  /  ALL0W(7) ,A|  ONGX < 7 ) . NPL OT . ZFPO . XM AX , TF I X  P0040 

DATA  (XT1TLF< T ) » 1=1. 1 0) /3«H  TIME  (RFC.)  ROORn 

DATA  (YTITLF(T).I  =  l.in)/3AH  PONOLINF  TFmPERATIJRF  (P>  P006O 

Al  ()NGV(  1  )=0.0  R0070 

no  io  1=1.7  Ronnn 

t  t  =  i  ponqo 

t  F  (  YL  T  M  -ALLOW  (  I  )  )  ?i),2n,10  ROlon 

10  CONTIMUF  ROlin 

wPITE  (6.1000)  YL IM  R012n 

1000  poHMATt///  37M  PPLOT  Will  NOT  BF  DONE  BFCAURF  YL T Mr  F12.R  ////)  P013O 

pF  TURN  ROUn 

20  YMAX  =AIlOW(IT)  P0150 

CALL  RSTFRM  R0160 

CALL  ORIt'ON(l?3. 1023, 24,024. IP. 1A. 5.5  )  R0170 

CALL  PLOtl  (1,1,ZFR0,XMAX,ZFR0»YmAX.X»Y»  NPLOT.l,  1H/  )  R01PO 

HO  30  1=1.6  P0190 

TAIL  I.APH  X  (ALONGX(t)  ,1  )  R0200 

CALL  l.ABM  Y  <  ALONPY  <  T  1  .  1  )  P0210 

30  aIONGY(I+1)  =  ALONGY(T)  +■  .2*  YMAX  P0220 

CALL  PRINT (2nn,P7R. 1 2. 0. 3P.XTTTLF)  R0230 

CALL  PRInT(47, 200.0. 12. 3P.YTITLF)  P0240 

CALL  PR TnT(123, 10 00.12.0. 72. TTTlF)  P02R0 

CALL  DMPMIF  P0260 

WFTUKN  R0270 

FMD  P0280 
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SIBFTC  CPIOT  F000O 

FIIRROIITTnF  CPI  OT  (X, Y1 » Y?.XLlM,Yl  TM1 ,YLTm2»TiTLF,  Y)  F001O 

nIMENFION  X(300) ,Yl(300) ,Y2(300) .YTITI  E( 1 0) »YY(?nOO) ,XTTTLE( 10)  F0020 

nlMFNFION  TITI  E(l?) *Y(30n) » ALONPy<7)  F0030 

UTMFNFIOn  CURVE(l) »VRUG<4) ,HPUG(7)  F0040 

COMMON  /  APC  /  ALL0W(7)*A|  ONPX < 7 ) , NPLOT » ?FPO , XmAX , TF I X  FOObO 

DATA  (VBUG(I) ,1=1,4)  /  100. 0,b0. 0,20.0, In. 0  /  F0060 

DATA  (HRUG(I)  ,1  =  1,7)  /  1  .  n  ,2 . 0 ,  5 .  0 * 1 0 . 0  ,  ?0 . 0  ,  SO .  0  r  100 . 0  /  F0070 

DATA  (  XT  1TLF  (  T  )  ,  1  =  1 , 1  0)  /3AH  TIME  (RFC.)  /  GOORO 

(lATA  (YTiTLF(T),I=l,tO)/^ftH  flTSTAMCE  (IN.)  /  FOOqn 

PATA  ONF/AHlnAO  /,TWo/4Hl460  /  FOlOO 

DATA  W0N/1H1  /,T00/1h?  /  F0110 

c  ***  four  (4)  CHArsCTF.PS  are  ALLOWED  FOR  CURVF(l)  F0120 

CURVE ( 1 ) =ONF  FO 130 

HFAC1R=HBUG( IFIX)  FO 1 40 

fympoi  =won  foiso 

YPIG  =AMay1  ( VLIM1  » Y|.TM?  )  F0160 

NCURVF  =1  F0170 

no  r  t=i »nplot  foiro 

1  Y Y ( 1 ) -  Yl(I)  F0190 

no  7  1=1,4  fo?oo 

Tl=  I  F0? 1 0 

IF  (  YBTG*100  .  -Al.LOWd)  )6,6,7  F0??0 

7  CONTINUE  F0?.3O 

WRITE  (6,1  OOn )  Yi.IMl.YLlM?  F0?40 

1000  FORMAT  (///  (PLOr  WTl  I  NOT  PF  DONE  RFCAUFF  Y| _IM1=F 1 ? . F , 1  pH  OR  F02FO 

1y(IM2=  F 1 ? , 6  ////  )  F0260 

RETURN  F0?70 

6  YNAX  =  ALLOW  (TD/100.  F0280 

VFAC I R=VbUG ( IT)  F0290 

CALL  RSTFPM  F0300 

CALL  OR TOON  (123, 102.3, 24.024, 18, 1R*S, 5)  F0310 

.1=1  F03?O 

70  no  10  I=J»NPLOT  F0330 

11=  I  F0340 

I Ml  =T-1  CP3S0 

IF(  YY(I)-Yd)  )20*ln»10  F03AO 

10  C.ONTINUF  F037O 

MOPT=NPLOT-J+1  F038P 

ll  =J  +  NoPT/2  F0390 

1VL0C=(Y|WAX-YV(LL) )*1«.*vFACTR  +24.  -4.  F0400 

1  ML  0C=  X(l.L)*l8.  /HF  aCTR  +123.  -48.  Cp410 

CALL  PRINT! INI  OC, IVLOC,  8 , 0 , 4  ,  Cl  IRVE )  F0490 

CALL  RL0T1  (1 ,1  ,?FRO, XMAX,7EP0,YMAX,X( J) ,YY( J) ,NOPT  ,1,SYmP0I)  F043O 

TF  (NCURVF-1  1  9(1, 85,90  F0440 

Pb  no  B6  I  =  l»NP|.oT  F04FO 

P6  YY ( I ) =Y2 ( T  )  F04 AO 

CURVE(1)=TW0  F0470 

FYMROLrTOO  FOAftn 

NCURVF  =  ?  F049O 

J=1  FOFOO 

GO  TO  70  F0F10 

20  NPT=IT-J  F0R2O 

I  L  =  J  +  ImPT /?  FOB30 

I VL0C= ( YM A  X— Y Y (LL) )  *  1  P . * Vc  ACTR  +24.  -4.  F0F40 

THLOC=  X ( I  L ) * 1 8 .  /HFaCTP  +123.  -48.  F0F50 

CALI  PR I NT ( I Hi  OC , T VLOC ,  A, 0,4, CURVE)  FOF60 


CAL  L  PL011  (1.1 .7FR0,xMAX,7ER0,YMaX,X(  J)  ,YY< J)  .NPT.l.RYMPOL)  50570 

no  50  IJ=  II,  NPLOT  50580 

JJ=  IJ  50590 

TP(YYdJ)-  Y(TJ)  150.ao.40  50500 

50  CONTINUF  50510 

I F  ( NCI  IRVt— 1 )  OO, 85.9,1  505?0 

40  J=  JJ  50550 

00  TO  70  50540 

90  fll  ON(iY ( 1  1=0.0  50550 

no  10O  1=1.5  50550 

CALI  LARKt X(AI  ONgX(I) ,  1)  50570 

CALL  LAPEI Y(AI  OMGY(T) ,1)  50580 

100  Al  ONGY<  1  +  1  >=AI  OA'GY(  T  )  +  ,?*YMAX  50590 

CALL  PR  I NT(2|I  0,975. 1p. 0,58, XTTTLF)  50700 

CALL  PRINT! 47, 200,0.  1 2 . 58 , YT I TLF 1  50710 

CALL  PRINT! 125, 1 000, 12, 0,72. TTTLF)  507?O 

CALL  riMPhllF  50750 

MFTUKN  50740 

FA'D  50750 


APPENDIX  D 


FORTRAN 

A 

AB 

ABLC 

ABLK 

B 

BB 

BL 

BLTEM 

BTEST 

C 

CB 

CHARC 

CHARK 

CP 

CPB 

CPC 

CPV 

CPX 

D 

DB 


PROGRAM  TERMINOLOGY 

Description 

"A"  coefficient  in  matrix,  single  subscript 

"A"  coefficient  in  matrix,  double  subscript 

specific  heat  of  material  at  TABL 

thermal  conductivity  of  material  at  TABL 

"B"  coefficient  in  matrix,  single  subscript 

"B"  coefficient  in  matrix,  double  subscript 

Total  thickness  of  backup  structure 

value  of  1460  isotherm  depth  from  previous  time  step 

test  to  determine  mode  of  heat  transfer  out  of  back  surface  of 
backup  materials 

"C"  coefficient  in  matrix,  single  subscript 

"C"  coefficient  in  matrix,  double  subscript 

specific  heat  of  material  at  TCHAR 

thermal  conductivity  of  material  at  TCHAR 

specific  heat  of  a  node  in  ablation  material 

specific  heat  of  backup  material  node 

specific  heat  values  in  char  specific  heat  table 

specific  heat  values  in  virgin  specific  heat  table 

specific  heat  values  in  backup  material  specific  heat  tables 

"D"  coefficient  in  matrix,  single  subscript 

"D"  coefficient  in  matrix,  double  subscript 

time  step  in  the  time  step  table 


DELTT 


FORTRAN 


Description 


DMP 

DRH0 

DT 

DTS 

DX 

DXB 

DXV 

DXX 

EMBB 

EMC 

EMFB 

EMV 

EMX 

END 

ERR1 

ERR2 

ERR3 

ERR4  _ 

FBL0W 

FC0NV 

FENV 

FRAD 

FTEST 


test  used  for  dumping  (DMP  =  0  skip  dump,  DMP  =  1.0  start  dumping) 
local  mass  flow  rate  of  ablation  gas 
time  step  from  the  time  step  table  in  hours 
time  step  from  time  step  table  in  seconds 
thickness  of  a  node  in  the  ablation  material 
thickness  of  a  node  in  a  backup  structure  material 

variable  ablation  node  thickness 

fixed  ablation  material  node  thi 

emissivity  of  back  surface  of  each  material  in  backup 
char  material  emissivity 

emissivity  of  front  surface  of  each  material  in  backup 
virgin  material  emissivity 

emissivity  of  front  surface  of  ablation  material 
code  word  for  plot  routine 

Control  numbers  for  printing  error  statements  when  an  input  or 
calculational  mistake  is  made 

blowing  efficiency  in  reducing  convective  heating 

factor  to  correct  convective  heating  rate  for  various  body  locations 
emissivity  -  view  factor  product  to  cabin  interior 

factor  to  correct  radiative  heating  rate  for  various  body  locations 

test  to  determine  mode  of  heat  transfer  into  front  surface  of 
backup  materials 
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FORTRAN 


Description 


FV 

G 

GAPX 

H 

H300 

HEAD 

HEADNG 

HEFV 

HTX 

HV 

HW 

HX 

IEM 

IPRC 

IPRCT 

IR 

IR1 

IR2 

NCASE 

NCPB 

NCPC 

NCPV 

NKC 


view  factor  for  external  environment 
defined  by  FORTRAN  statement 
gap  width  between  backup  materials 
film  coefficient  between  backup  materials 
enthalpy  of  air  at  300°  K 

any  72  alphanumeric  characters  used  to  identify  problems  being 
run  -  printed  at  top  of  first  page  of  output 

any  72  alphanumeric  characters  used  to  identify  each  input  section 
film  coefficient  to  cabin  environment 
total  enthalpy 

heat  of  degradation  of  virgin  material 

wall  enthalpy  computed  from  enthalpy  —  temperature  table 
enthalpy  values  in  enthalpy  table 

test  used  to  determine  if  front  surface  is  virgin  or  char  for  using 
proper  emissivity 

variable  print  frequency  in  time-step  table 
present  print  control  number 

test  to  determine  if  node  temperature  is  greater  than  TABL 
test  used  in  determining  node  density  at  TX1  temperature 
test  used  in  determining  node  density  at  TX2  temperature 
number  of  problems  to  be  run 

number  of  points  in  each  backup  material  specific  heat  table 
number  of  points  in  char  specific  heat  temperature  table 
number  of  points  in  virgin  specific  heat  temperature  table 
number  of  points  in  char  thermal  conductivity  -  temperature  table 


FORTRAN 


Description 


NKPB 

RKV 

NMB 

NP 

NPBS 

NPF 

NPL0T 

NPM 

NHP 

NPTT 

NREC 

NTIIAPT 

NXA 

MB 

MC  \ 

NXD 

ME 

y 

QBL0CK 

qcjzSn 

QC0M 

QHW 

QIN 

QL0SS 

Q0XID 


number  of  points  in  each  backup  material  thermal  conductivity  table 

number  of  points  in  virgin  thermal  conductivity  temperature  table 

number  of  materials  in  backup  structure 

number  of  node  points  in  ablation  material 

total  number  of  node  points  in  backup  structure 

total  number  of  points  in  heat  shield  structure  (NP  +  NPBS) 

output  plot  control  number 

number  of  nodes  per  material  in  backup 

number  of  points  in  enthalpy  -  temperature  table 

number  of  points  in  time-step  table 

number  of  points  in  surface  recession  —  temperature  or  time  table 
number  of  points  in  trajectory  input  table 


dummy  indices  for  subroutine  SAVE 


amount  of  convective  heat  blocked  due  to  mass  injection  into 
boundary  layer 

trajectory  table  convective  heating  rates 
cold  wall  convective  heat  rate  at  present  time  step 
hot  wall  convective  heat  rate  without  blowing 
net  heat  flux  into  front  surface 

boundary  condition  for  heat  transfer  to  cabin  interior 
heating  rate  due  to  combustion 
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FORTRAN 


Description 


QRAD 

QRADX 

QUIT 

R 

RBI 

RB2 

RH0 

RH0BX 

RH0C 

RH0V 

RH0Y1 

RH0Y2 

S 

SD0T 

SAVEIT 

SAVEXX 

SAVX 

SAVY1 

SAVY2 

SAVY3 

SAVY4 

SAVY1X 

SAVY2X 


trajectory  table  radiative  heating  rates 
radiative  heat  flux  at  present  time  step 
code  word  for  plot  routine 

thermal  resistance  due  to  conductivity  between  nodes  in  the  ablation 
material 

thermal  resistance  due  to  conductivity  between  past  and  present  node 
in  backup  material 

thermal  resistance  due  to  conductivity  between  present  and  forward 
node  in  backup  material 

density  of  an  ablation  material  node 

density  of  individual  materials  in  backup 

mature  char  material  density 

virgin  ablation  material  density 

density  of  node  at  past  time  step 

density  of  node  at  present  time  step 

thermal  capacity  of  a  node  in  the  ablation  material 

surface  recession  rate 

depth  of  1060  isotherm  at  any  given  time 
time  corresponding  to  maximum  depth  of  1460  isotherm 
time  corresponding  to  maximum  depth  of  1060  isotherm 
surface  recession  depth  at  maximum  1060  isotherm  depth 
bondline  temperature  at  maximum  1060  isotherm  depth 
term  that  will  contain  maximum  depth  of  1060  isotherm 
depth  of  1460  isotherm  at  maximum  1060  isotherm  depth 
surface  recession  depth  at  maximum  1460  isotherm  depth 
bondline  temperature  at  maximum  1460  isotherm  depth 
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FORTRAN 


SAVY3X 

SAVY4X 

SR 

T 

TABL 

TCHAR 

TCP 

TCPC 

TCPV 

TEMDI 

TEMPI 

TENV 

TEST2 

TDMP 

TIME 

TINT 

TITLE 

TKC 

TKV 

TL 

TLIM 

TREC 

TS 

TTABLE 


Description 

depth  of  1060  isotherm  at  maximum  1460  isotherm  depth 
term  that  will  contain  maximum  depth  of  1460  isotherm 
surface  recession  values  in  surface  recession  table 
present  time 

temperature  at  which  ablation  starts 
temperature  at  which  ablation  stops 

temperature  values  in  backup  material  specific  heat  tables 
temperature  values  in  char  specific  heat  table 
temperature  values  in  virgin  specific  heat  table 
arbitrary  initial  temperature  distribution  values 
constant  initial  temperature  distribution  value 
interior  cabin  temperature 

test  to  determine  proper  heat  shield  initial  temperature  distribution 

time  to  start  dumping  or  printing  information  used  in  checkout  of 
program  (sets  DMP  =  1.0) 

trajectory  table  time  values 

starting  time  of  problem 

control  card  used  for  reading  in  new  data  for  successive  problems 
temperature  values  in  char  thermal  conductivity  table 
temperature  values  in  virgin  thermal  conductivity  table 
total  thickness  of  heat  shield  structure  (VL  +  BL) 
time  limit  of  problem 

surface  temperature  or  time  at  which  char  removal  is  to  start 
temperature  or  time  values  in  surface  recession  table 
time  values  in  time -step  table 
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FORTRAN 


Description 


TTUL  equals  TUL  if  VPT  =  0  or  equals  TX2  if  VPT  =  1  -  used  in  computing 
char  properties 

TUL  maximum  value  of  TX1  and  TX2 

TUL1  maximum  TX1  values  -  used  in  computing  gas  ablation  rate 

TUL2  maximum  TX2  values  -  used  in  computing  gas  ablation  rate 

TV  sink  temperature  of  external  environment 

TW  temperature  values  in  enthalpy  table 

TX1  temperature  of  nodes  at  past  time  step 

TX2  temperature  of  nodes  at  present  time  step 

TX2C  temperature  at  fixed  locations  in  ablation  material  as  defined  by  XC 

TX2T  temporary  storage  of  TX2  temperatures  for  computing  thermal  properties 

TXK  temperature  values  in  backup  material  thermal  conductivity  tables 

TX0  initial  temperature  at  front  surface  of  heat  shield  for  computing 

linear  temperature  gradient 

TY  temperature  distribution  at  forward  time  step 

VEL  trajectory  table  velocity  values 

VELX  trajectory  velocity  at  present  time  step 

VL  initial  virgin  material  thickness 

VLI  initial  ablation  material  thickness 

VLTEM  value  of  1060  isotherm  depth  from  previous  time  step 

VLV  variable  ablation  material  thickness 

VPT  test  to  determine  if  properties  are  irreversible  with  temperature 

WEKEEP  depth  of  1460  isotherm  at  any  time 

XBM  thickness  of  individual  materials  in  backup 

XC  fixed  location  of  nodes  in  the  ablation  material 

XI  node  number 
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FORTRAN 


Description 


XIDNT 

XK 

XKB 

XKC 

XKV 

XL0ST 

XLSTI 

XLSTV 

XMDC 

XMDG 

XMD0 

XMDT 

XNP 

XNPM 

XPL0T 

XV 

YK 

YPL0T1 

YPL0T2 

YPL0T3 

YPL0T4 

ZZZ 


any  72  alphanumeric  characters  to  identify  each  material 

thermal  conductivity  values  in  backup  material  thermal  conductivity 
table 

thermal  conductivity  of  backup  material  node 

thermal  conductivity  in  char  thermal  conductivity  table 

thermal  conductivity  value  in  virgin  thermal  conductivity  table 

amount  of  solid  ablation  material  lost  in  a  time  step  due  to  surface 
movement 

distance  from  original  surface  to  present  front  surface  location, 
inches 

distance  from  original  surface  to  present  front  surface  location, 
feet 

mass  loss  rate  of  char 

mass  gas  ablation  rate  due  to  pyrolysis  of  virgin  material 

mass  flux  rate  of  oxygen  to  surface 

total  ablation  rate 

number  of  nodes  in  ablation  material 

number  of  nodes  per  backup  material 

time  to  be  written  on  tape  and  plotted 

location  of  nodes  in  variable  ablation  material  thickness 

thermal  conductivity  of  a  node  in  ablation  material 

recession  depth  to  be  written  on  tape  and  plotted 

bondline  temperature  to  be  written  on  tape  and  plotted 

1060  isotherm  depth  to  be  written  on  tape  and  plotted 

1460  isotherm  depth  to  be  written  on  tape  and  plotted 

ratio  to  determine  when  the  limiting  value  of  heat  blockage  has  been 
reached 
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TABLE  I.-  SAMPLE  PROBLEM  IRFOT 


NOTE;  WRITE  NUMBERS  10  ,  LETTERS  I  <f>  U&H  G  SYts/tBOi-S  /■ 


TABLE  I.-  SAMPLE  PROBLEM  INPUT  -  Concluded 


note:  WRITE  NUMBERS  10  ,  LETTERS  I <f> U&ZC  j  SYMBOLS  /•■>  * 


TABLE  I.-  SAMPLE  PROBLEM  INPUT 


(t>)  Fortran  data  card  listing 


1 

TYPICAI  OARRING  ABLATOR  -  TEST  pare  -  4/6/65  DONALD  M.  ClIPPY 


+bOO*U 

+  1)0 

+  11.0 

hid 

? 

1 

til,  u 

+  U" 

+  11 , 1 

+  110 

1  oo 

+  b  0  0  •  U 

+  00 

+  0 . 1 

+  00 

1  00 

+  1.0 

+  00 

+  1.0 

+  00 

TYPICAL  1 

CH».RKInG 

A  HI  AT  I  ON  MATFRTAI.  PROPFRTIFS 

+  00 

+  lAbll.O 

+  0  0 

+  U.  0 

+oo 

+  34 . 0 

+  U0 

+?o.o 

+  00 

+  0.00 

+  00 

+  L1  .65 

+  11  o 

+  0.7S  • 

+  0  0 

+1?9.06 

+00 

+  1.50 

+  00 

+250.0 

+  0  0 

+  0.0 

+00 

+  1  .U 

+  0  0 

TI..0 

+  00 

+  0.12 

+oo 

+  0.43 

+  00 

+0.070 

+  00 

+  0.43 

+00 

31 

2 

2  ° 

2  2 

+  i<i60 .  n 

+u  o 

+  0.1? 

+  00 

+  1.0 

+  U4 

+  0.1? 

+  00 

+1460.0 

+  00 

+  0.43 

+  00 

+  1.0 

+  04 

+  0.43 

+  00 

+  360  » 1.1 

+  00 

+  o  .  06s 

+  00 

+400.0 

+  00 

+0 . 065 

+  00 

+560.0 

+  00 

+0.0656 

♦  On 

+to6U.U 

+  oO 

+0 . 06t> 

+  00 

+  /AO.  0 

+  0  0 

+0.067? 

+  00 

+R60.0 

+  00 

+0.0684 

+  00 

+960.0 

+  0  o 

+  o .  U69 

+  00 

+1060.0 

+  00 

+0.070 

+  00 

+1160.0 

+  00 

+0.070 

+  00 

+  360 . u0 

+  (|0 

+  0.43 

+  00 

+1100,0 

+  00 

+  0.43 

+  00 

+  0.0 

+  00 

+y .  o 

-04 

+  fc>0  o .  0 

+  00 

+  9.0 

-04 

NO  TRAJECTORY  -  0=95 

RTU/Sfc  C-SQFT 

£- 

+  0.0 

+  l.l  0 

+ys.  n 

+  U0 

+  0.0 

+  00 

+2 . 9?5 

+  04 

+600. U 

+u  o 

+  95,0 

+  00 

+  .I.0 

+00 

+2 . 9?5 

+  04 

1 

3 

+  0,1 

+  00 

+  3.0 

+  11  0 

9 

2 

HACK  UP  MAT t  R I Al 

0.1 

L  INCHES 

THTCK 

+360. U 

+  0  o 

+ti  »  ribs 

+  1)0 

+460 • 0 

+00 

+0.065 

+  00 

+560.0 

+  00 

+0.06R6 

+  00 

+  r6U  •  0 

+  00 

+  l» • 0bfo 

+  00 

+  76  u .  0 

+  00 

+0 . U67? 

+  00 

+860.0 

+  00 

+0.0684 

+  00 

+960.0 

+  00 

+  1-.069 

+  0  0 

+1060.0 

+  00 

+  0.07 

+  00 

+1160.0 

♦  oo 

+  0,07 

+  00 

+360 .00 

+  1)0 

+  U  ,43 

+  00 

+1100.0 

+  00 

+  0.43 

+  00 

+  34.0 

+  00 

+  U  •  1 

+  iio 

+  0.9 

+oo 

+  0.9 

+  00 

+  U  .  0 

+  0" 

.  .1 ,  U 

+  oo 

+  0.0 

+  00 

+  0.0 

+  00 

heat  TkAnSFFk  TO  PAPIN 

F'\|VIP(>MMI-N  1 

-  HENvrO.O 

+660,0 

+  U0 

+  0,0 

+  1)0 

+  0.0 

+  00 

+  0.0 

+  00 

INITIAL  TPMP+RATuPE  IS 

CONSTANT 

+  0.0 

+  II  n 

+63G.U 

+o0 

+  ">30.0 

+  0  0 

U2 

+  U  .  0 

+  nn 

+i. .  n 

+  0  0 

+342.9 

+00 

+1400.0 

+  00 

+449.7 

+  no 

+1800.0 

+  00 

+  ol  7.2 

+  oo 

+  2400. n 

+  00 

+791.0 

+  00 

+3000.0 

+  00 

+978.0 

+  00 

+3600 . 0 

+  00 

+1113.0 

+  11 0 

+4000.0 

+  00 

+1P0O.O 

+  00 

+4924 . n 

+  on 

+1 300.0 

+  00 

+4486.0 

+  00 

+  1 A  0  O , 0 

+  0  0 

+  4  '23.  n 

+  0  0 

+1600.0 

+  00 

+4936.0 

+  00 

+1600.0 

♦  00 

+5127.0 

+  00 

+1700.0 

+  0  0 

+s?yo. o 

+  00 

+1000 .0 

+  00 

+  54  54.0 

+oo 

+1°00. 0 

+  1)0 

+5696.0 

+  00 

+2000.o 

+  0" 

+6728.0 

+  00 

+2100.0 

+  00 

+5851.0 

+  00 

+2?on.n 

+0O 

+5068.0 

+  00 

+  2300.0 

+  00 

+6078 . 0 

+  0  0 

+2400.0 

+nn 

+6106.0 

+  00 

+2600.0 

+oo 

+6?91 .0 

+  00 

-.2P00 . 0 

+on 

+6395.0 

+  00 

+2700.'> 

+  00 

+6497.0 

+  00 

+2800. 0 

+  00 

+6597.0 

+  00 

+2°00,0 

.+  1)0 

+6699.0 

+  0  0 

+3000. 0 

+00 

+6805.0 

+  00 

+3100.0 

+  0O 

+6918.0 

+00 

TABLE  I  -  SAMPLE  PROBLEM  INPUT  -  Concluded 


(b)  Fortran  data  card  listing 


+1200. 
+35O0. 
+  3hOo  . 
+4  mo. 

+  U4  n  o • 


+  Un  +  7i|C'U  .  u 
+  (ifi  +  7uP0 ,  0 
+U"+/470.0 
+ijo+r^oo .  ii 
+uo+Q(.nu.ii 


+UO+3300.  n 
+00+3500.0 

+itn+.i4n(i.ii 

+uo+42on. u 
+oo+ubno.n 


+  1)0+7175.0 
+nn+7bio.o 
+on+pi?o.o 
♦no+A700.o 

+  00+9150. o 


+0O+14O0.0 
+00+1700.0 
+  00+40  0  0.0 
+00+4300.0 
+00+4b00.0 


+00+7350.0  +00 
+0O+7BO0.0  +00 
+00+0300 . 0  +00 
+00+BBR0.0  +00 
+00+0270 . 00  +00 


TABLE  II.-  SAMPLE  PROBLEM  OUTPUT 


TYPICAL  CHARRING  ARI  ATOR  -  TEST  CASE  -  4/6/65  DONALD  M.  CURRY 


INPUT  DATA. 


TIME  LIMIT=6.D000E  0?  INITIAL  TIME=0. 


NPTT=  2 


TIME 

0. 

6.0000E  02 


TIME  STEP 
l .OOOOE-Ol 
t.OCOOE-OI 


PRINT  CONTROL 
100 
100 


FCON V=  l.OOOOOE  00  FRAI)=  l.OOOOOE  00 


TYPICAL  CHARRING  ABLATION  MATERIAL  PROPERTIES 


T  ABL=  1 .06000 E 

03  TCHAR=  1.46000b  03  TREC= 

0. 

RHOV=  3.40000E 

01 

RHOC  = 

2.00000E 

01 

FBLOW-  0. 

EMV=  6 • 50000E-0 1  EMC= 

7 • 50000E-0 l 

N300=  1.29060E 

02 

VL  = 

1.50000E 

00 

H V=  2.6 

UOOOD 

02  VPT=  0.  F  V  = 

1.000006  00 

T  V=  0. 

CHARK- 

1.20000E- 

01 

C  HAR  C=  4 . 30000E 

-01  ABLK=  7.00000E-C2  ABLC= 

4 • 30000E-01 

NP=  31 

MKC 

=  2  N  C  P  C  =  2  NKV=  9 

NCPV=  2 

NREC=  2 

VIRGIN  MATERIAL 

THERMAL 

SPECIFIC 

TEMPERATURE 

CONDUCT  I VI TY 

TEMPERATURE 

HEAT 

3.60000E 

02 

6  •  60000E-02 

3.60000E  02 

*♦ . 30000  E-0 1 

4.60000E 

02 

6. 50000E-02 

1.10000E  03 

4  •  30000E-0 1 

5.60000E 

02 

6 . 55000E-U2 

6. 6 0000 E 

0? 

6 • 60000E-02 

7.60000C 

02 

6.7200OE-02 

8.60000E 

0? 

6.04000E-02 

9.60000E 

02 

6.90000E-02 

I • 06000b 

03 

7.000006-02 

1.L6000E 

03 

7 • OOOOuE-02 

CHAR  MATERIAL 

THERMAL 

SPECIFIC 

TEMPERATURE 

CONDUCTIVITY 

TEMPERATURE 

HEAT 

1 • 46000  E 

0  3 

l • 20 009 £-01 

1.46000E  03 

•t  .  JOUOut  -VII 

l.OOOOOE 

04 

1.20000E-01 

l.OOOOOE  04 

4 . 30000E-0 1 

SURFACE  RECESSION  TABLE 


TIME 

0. 

6.00000E  02 


SK  -  1 N/ SEC 
9. 0C000E-04 
9 . 0C000E-04 


NO  TRAJECTORY  -  0=95  BTU/SEC-SQFT 


NO.  OF  TRAJECTORY  POINTS 


2 


TABLE  II.-  SAMPLE  PROBLEM  OUTPUT  -  Continued 


TIME 

0  CONVECTIVE 

Q  RADIATIVE 

VELOCITY 

0. 

V.t>UUU0t  Ui 

U. 

2. v^bout  i 

6 . OOOOOE  02 

9.50000F  01 

0. 

2.92500E 

PROPERTIES  OF  BACKUP  STRUCTURE 


NO.  CF  MATERIALS  IN  BACK-UP  SHIELD=  1 
TOTAL  MUMRER  OF  NUDES  IN  BACK-UP  SH I  ELD  =  3 

THICKNESS  OF  BACK-UP  SHI  EL 0=  l.OOOOOE-Ol 


BACKUP  MATERIAL 


TMP^MAI 


TEMPERATURE 

CUNUUL l 1 V i  IT 

3.6000'JE 

02 

6. 5 OOOOE-O 2 

4.6000CF 

02 

6. 5C00 05-02 

5. 60000 E 

02 

6. S 5000 E- 02 

6 • 60000E 

02 

6 . 6C0UOE-02 

7. 6  00  OOF. 

o; 

6. 7 2000c- 02 

8. 60000 E 

02 

6.  a  40  0011-02 

9. 600 00 E 

02 

6.9000-JE-0? 

l. 06000 E 

03 

7. OOOOOE- 02 

l.  I6000C 

03 

7 • 00000E-02 

MATERIAL 

1 

DENSITY 

3.4000E  01 

0.1  INCHES  THICK 


TEMPERATURE 
3.60000E  02 
1.10000E  03 


EMISSI 

THICKNESS  FRONT 

1 • OOOOE-OI  9 • 0000c-01 


SPECIFIC 
nE  AT 

4 • 30000E-0 1 
4 • 30000E-0 1 


VITY 

BACK  NODES/MATERIAL 

9 . OOOOE-O 1  3.0000E  00 


ADDITIONAL  DATA  FUR  INDIVIDUAL  MAlERlALS  IN  BACKUP  STRUCTURE 


MATERIAL 

I 


FILM  COEFFICIENT 

0. 


GAP  THICKNESS 

0. 


FTEST 

0. 


BTEST 

0. 


HEAT  TkANSFER  TO  CABIN  ENVIRONMENT  -  HENV=0.0 
T  EMPtK  ATURE=  b  •  60000E  02  FILM  COEFFICIENT  C.  VIEW  FACTOR=  0. 


INITIAL  TEMPERATURE  IS  CONSTANT 

TEMPFRATURF  DISTRIBUTION  IN  HEAT  SHIELD  IS  UNIFORM  AND  EQUAL  TO  5.3000F  02 


Q  LO S T=  0. 
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TABLE  II.-  SAMPLE  PROBLEM  OUTHJT  -  Concluded 


OUTPUT  DATA. 


TIME-  9.90000E  00  OCUNVECTIVE*  9.50000E  01  QRAOIATIVE*  0.  VELOCITY*  2.92500E  04 

GAS  ABLATION  RATE*  0.  CHAR  ABLATION  R ATE=  5.40000E  00  TOTAL  ABLATION  RATE=  5.40000E  DO 

RECESSION  DEPTH*  9.00000E-03  OHOT  WALL*  8.99282E  01 

TEMPERATURE  CIS  TR IbUT I  UN  IN  HEAT  SHIELD  AT  THE  END  OF  THE  TIME  STEP,  T=  l.OOOOOE  01  SECONDS 


TEMPERATURE  DISTRIBUTION  IN  THE  ABLATING  MATERIAL 


3.7 902 2E  03  2.33046E  03  1.06300E  03  6.38075E  02  5.47600E  02  5.32434E  02 
5.30289E  02  5.30030E  02  5.30002E  02  5.30000E  02  5.30000E  02  5.30000E  02 
5 • 30000E  02  5.30000E  02  5.29999E  02  5.30000E  02  5.30000E  02  5.3DOOOE  02 
5. 300006  02  5.30000E  02  5.30000E  02  5.30000E  02  5.30000E  02  5.30000E  02 
5.30000E  02  5.30000E  02  5.30000E  02  5.29999E  02  5.30000E  02  5.30000E  02 
5 #299996  02 


TEMPERATURE  DISTRIBUTION  IN  THE  BACK-UP  STRUCTURE 


5.29999E  02  5.30000E  02  5.30000E  02 


TIME*  1 • 990006  01  QCONVECTIVE*  9.50000E  01  QRAOIATIVE*  0.  VELOCITY*  2.92500E  04 

GAS  ABLATION  RATE*  2.90720E  01  CHAR  ABLATION  RATE*  5.40000E  00  TOTAL  ABLATION  RATE*  3.44720F  01 
RECESSION  DEPTH*  1.80000E-02  QHOT  WALL*  8.99173E  01 

TEMPERATURE  DISTRIBUTION  IN  HEAT  SHIELO  AT  THE  END  OF  THE  TIME  STEP,  T=  2.00000E  01  SECONOS 


TEMPERATURE  DISTRIBUTION  IN  THE  ABLATING  MATERIAL 


3.80935E  03  2.87678E  03  1.645756  03  9.77202E  02  A.72698E  02  5.66204E  02 
5.37857E  02  5.  31496E  02  5.30254E  02  5.30038E  02  5.30005E  02  5.3DD006  02 
5.29999E  02  5.299996  02  5.299996  02  5.29999E  02  5.29999E  02  5.29999E  02 
5.29999E  02  5.29999E  02  5.29999E  02  5.29999E  02  5.29999E  02  5.29999E  02 
5.29999E  02  5.29999E  02  5.29999E  02  5.29999E  02  5.29999E  02  5.29999E  02 
5 • 29999E  02 


TEMPERATURE  DISTRIBUTION  IN  THE  BACK-UP  STRUCTURE 


5.29999E  02  5.29999E  02  5.29999E  02 


TIME*  2.99000E  01  QCONVECTIVE*  9. 50000E  01  QR ADIAT IVE=  0.  VELOCITY*  2. 92500E  04 

GAS  ABLATION  RATE*  1.253306  01  CHAR  ABLATION  RATE*  5.40000E  00  TOTAL  ABLATION  RATE*  1.79330E  01 
RECESSION  DEPTH*  2. 700006-02  OHOT  WALL*  8.98427E  01 
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Figure  2.  -  Schematic  diagram  of  charring  ablator  thermal  protection  system. 
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Figure  9.  -  Plot  program  bond  line  temperature  curve  from 
typical  charring  ablator  test  case. 
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Figure  10.  -Plot  program  1060°R  and  146 0°R  isotherm  curves  from 
typical  charring  ablator  test  case. 
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